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Until  recently,  the  classical  Streeter-P helps  equation  was 
widely  used  in  predicting  the  dissolved  oxygen  deficit  to  organic 
waste  loads  discharged  into  a  stream.  However,  the  wide  deviation 
of  predicted  values  of  oxygen  deficit  from  the  observed  value  lead 
several  investigators  to  modify  the  S tree ter -Phelps  formulation  by 
taking  into  account  other  factors  like  benthal  demand,  photosynthesis, 
algal  respiration,  etc.,  which  affect  the  dissolved  oxygen  in  the 
stream,  in  order  to  narrow  the  gap  between  the  observed  and  the  pre- 
dicted values.  Other  investigators  postulated  statistical  models 
for  predicting  critical  oxygen  concentrations  using  stream  flow, 
temperature,  and  5-day  BOD  as  independent  variables,  mainly  to  avoid 
the  errors  involved  in  the  predictions  due  to  the  fluctuations  in 
the  values  of  deoxygenation  and  reaeration  coefficients. 

This  work  attempts  to  predict  dissolved  oxygen  deficits  in  a 
stream  with  known  initial  conditions  by  taking  into  account  the  varia- 
tions in  deoxygenation  and  reaeration  coefficients.  A  hypothetical 
stream  situation  is  used  to  establish  the  significance  in  predicting 
dissolved  oxygen  deficit.  Statistical  models  are  formulated  and 
tested  for  the  variations  in  these  coefficients  using  published  data. 
Simulation  techniques  using  the  Monte  Carlo  method  are  employed  in 
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predicting  the  probabilistic  variation  in  dissolved  oxygen  deficits 
for  known  initial  conditions  and  the  results  are  verified  with  the 
survey  data  observed  for  the  Ohio  River-Cincinnati  Pool  reach.  The 
predicted  results  using  probabilistic  model  are  found  to  agree  with 
the  observed  values  within  practical  limits  and  give  more  consistent 
results  than  conventional  methods. 
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I.   INTRODUCTION 

The  rapid  expansion  of  population  and  industrialization 
has  resulted  in  increasingly  difficult  problems  of  water  resources 
management.  The  most  critical  of  these  problems  is  the  protection 
of  water  resources  from  the  ravages  of  pollution  by  the  discharge  of 
wastes  which  are  increasing  in  both  volume  and  complexity.  There  is  a 
general  awakening  and  a  demand  that  something  must  be  done  about  stream 
pollution.  The  daily  news  media  report  various  aspects  of  water  pol- 
lution with  considerable  frequency  and  in  the  past  decade  many  national 
laws  have  been  enacted  dealing  with  Federal  participation  in  pollution 
control  activities. 

If  a  pollution  control  program  is  to  be  both  effective  and 
economical,  engineers  must  possess  the  ability  to  predict  the  effects 
of  specific  waste  discharges  on  the  environment.  Without  this  know- 
ledge, administrative  determination  of  the  required  degree  of  treatment 
for  waste  discharges,  existing  or  proposed,  can  only  be  speculative. 
A  low  estimate  of  the  required  degree  of  treatment  will  result  in  un- 
desirable conditions,  while  an  overestimate  will  create  unjustifiable 
economic  burdens  on  the  waste  discharger  without  commensurate  benefit 
to  the  environment. 

The  determination  of  acceptable  levels  of  water  quality  is 
usually  the  task  of  a  regulatory  agency.  These  governmental  agencies 
have  used  several  approaches  to  control  the  water  quality  in  streams. 
One  approach  is  through  the  use  of  treatment  standards  which  require 
specified  degrees  of  treatment,  usually  in  terms  of  biochemical  oxygen 
demand  (BOD)  reduction,  suspended  solids  and  coliform  removal.  Effluent 


standards  have  also  been  established  wlich  limit  the  concentrations 
of  the  various  constituents  in  the  effluent  released  into  the  stream. 
These  regulatory  activities  are  attuned  to  the  concept  of  equity,  which 
require  the  same  degree  of  removal  of  BOD  and  solids  from  all  waste 
sources  irrespective  of  the  volume  of  waste  (Jacobs,  1965).  These 
standards  are  used  because  they  are  relatively  easier  to  administer, 
but  each  one  deals  only  indirectly  with  the  basic  problem  of  stream 
quality. 

The  establishment  of  stream  standards  requires  the  polluter 
to  regulate  his  effluent  such  that  at  least  a  minimum  level  of  stream 
quality  is  maintained  at  all  times*  Though  the  aim  of  all  the  three 
types  of  standards  is  to  assure  acceptable  levels  of  stream  quality,  it 
is  directly  accomplished  only  in  the  case  of  stream  standards.  The 
adoption  of  stream  standards  as  a  tool  for  regulatory  control  of  pollu- 
tion is  not  without  limitations.  In  the  first  place  it  is  very  difficult 
to  administer  and  enforce,  especially  if  there  is  more  than  one  polluter. 
Also,  the  application  of  rigid  stream  standards  to  large  areas  may  well 
become  a  barrier  to  orderly  economic  development,  thus  defeating  some 
of  the  benefits  to  be  derived  from  the  equitable  use  of  water  resources 
within  a  given  locality  or  region  (Jacobs,  1965). 

The  Important  point  which  needs  to  be  made  here  is  that  the 
concept  of  stream  standards  and  effluent  standards  is  not  mutually  ex- 
clusive. In  most  cases  both  are  necessary.  For  example,  the  Illinois 
Sanitary  Water  Board  (1966)  has  adopted  rules  and  regulations  pertaining 
to  sewage  and  industrial  waste  treatment  requirements,  effluent  criteria, 
and  water  quality  criteria  for  lakes  and  rivers,  in  order  to  protect  the 
water  resources  of  the  state.  It  is  true  that  the  amount  of  weight 


placed  on  a  particular  type  of  standard  depends  on  a  given  situation, 
the  type  of  waste  and  the  stage  of  development  of  the  pollution  abate- 
ment program.  As  pollution  abatement  progress  is  made  and  as  the  de- 
mands for  various  water  uses  grow,  more  and  more  emphasis  will  have  to 
be  given  to  the  stream  standards. 

Though  the  stream  quality  standards  specify  minimum  accept- 
able levels  for  a  variety  of  stream  quality  parameters  like  bacteria, 
dissolved  solids,  chemical  constituents,  etc.,  the  parameter  most  com- 
monly used  as  a  measure  of  the  pollution  from  biodegradable  waste,  by 
investigators  concerned  with  stream  sanitation  aspects,  is  the  stream's 
dissolved  oxygen  (DO)  concentration.  In  this  study,  the  primary  atten- 
tion will  be  directed  to  this  parameter,  though  it  is  conceptually 
feasible  to  extend  the  ideas  to  deal  with  other  parameters  as   well. 

With  the  enactment  of  the  Water  Quality  Act  of  1965,  state 

regulatory  agencies  either  have  adopted  or  are  in  the  process  of  adopting 

(as  of  May  1968)  standards  for  interstate  waters.  Among  other  parameters 

considered  for  stream  standards,  the  Illinois  Sanitary  Water  Board 

(1966)  has  adopted  different  standards  for  dissolved  oxygen  in  Illinois 

rivers,  depending  on  the  water  use.  Thus,  the  criteria  adopted  for 

the  aquatic  life  sector  of  streams  is: 

For  maintenance  of  well  balanced  fish  habitats  the 
dissolved  oxygen  content  shall  be  not  less  than  5.0 
mg/1  during  at  least  16  hours  of  any  24  hour  period, 
nor  less  than  3.0  mg/1  at  any  time, 

and  for  the  industrial  water  use  sector,  It  is: 

Not  less  than  3.0  mg/1  during  at  least  16  hours  of 
any  24  hour  period,  nor  less  than  2.0  mg/1  at  any 
time. 


Though  the  standards  adopted  for  dissolved  oxygen  are  defini- 
tive in  nature,  there  is  a  finite  probability  that  the  dissolved  oxygen 
is  likely  to  fall  below  the  set  standard  during  16  hours  of  any  24-hour 
period  and  the  standard  set  for  the  minimum  dissolved  oxygen  at  any 
time.  There  is  a  growing  realisation  among  the  investigators  for  the 
need  of  probabilistic  stream  standards  (loucks,  1967;  Ledbetter  and 
Gloyna,  1964;  Thayer,  1966).  Loucks  (1967)  proposes  a  probabilistic 
stream  standard  for  dissolved  oxygen  as  follows: 

The  dissolved  oxygen  concentration  in  the  stream  during 
any  7  consecutive  day  period  must  be  such  that 

1.  The  probability  of  its  being  less  than  4  mg/1 
for  any  1  day  is  less  than  0,2;  and 

2.  The  probability  of  its  being  less  than  2  mg/1 
for  any  1  day  is  less  than  0,1  and  for  any  2  or  more 
days  is  less  than  0.05. 

OBJECTIVE  OF  THE  STUDY 

In  order  to  design  any  type  of  treatment  plant,  it  is  necessary 
to  know  what  volume  and  concentration  of  pollutant  may  be  discharged  to 
the  stream  so  that  the  stream  standards  are  not  violated.  This  implies 
that  there  should  be  methods  to  predict  the  response  of  the  receiving 
stream  to  waste  loads  placed  in  it.  Several  mathematical  formulations 
are  available  for  predicting  the  dissolved  oxygen  responses,  most  of 
which  are  based  on  the  pioneering  work  of  Streeter  and  Phelps  (1925). 
In  all  these  formulations,  the  reaction  velocity  coefficients  affecting 
the  rates  of  BOD  removal,  atmospheric  reaeration,  etc.,  are  taken  as 
constants,  though  it  was  recognised  by  a  few  investigators  (LeBosquet 
and  Tsivoglou,  1950;  Eckenfelder  and  O'Connor,  1961)  that  they  are  far 
from  being  constants.  The  main  aim  of  this  work  is  to  ascertain  the 
significance  of  the  variations  in  the  velocity  coefficients  K}  and  K2 


in  defining  the  dissolved  oxygen  (DO)  response  of  a  receiving  stream 
and  to  develop  a  procedure  for  determining  the  DO  taking  the  variability 
in  these  rate  coefficients  into  consideration,  if  these  variations  are 
significant. 

The  specific  objectives  of  this  study  are: 

1.  To  examine  the  relative  importance  of  the  variations  in 
the  reaction  velocity  coefficients  affecting  the  bacterial  oxidation 
of  the  organic  matter  and  the  atmospheric  reaeration  of  river  water, 
in  predicting  the  dissolved  oxygen  responses  of  the  receiving  stream. 

2.  To  determine  the  nature  of  variations  of  these  velocity 
coefficients  and  to  formulate  and  test  the  hypotheses  concerning  their 
chance  variations. 

3.  To  develop  a  procedure  for  predicting  the  dissolved 
oxygen  in  a  river  downstream  of  a  waste  source  by  taking  into  account 
the  variations  in  these  velocity  coefficients. 

4.  To  enumerate  quantitatively  the  chance  variations  in  dis- 
solved oxygen  responses  in  streams  in  terms  of  probability  measure. 

The  last  of  the  four  objectives  mentioned  above  is  extremely 
significant  in  the  light  of  the  guidelines  established  for  water  quality 
standards  under  the  Water  Quality  Act  of  1965,  wherein  it  is  stated  that 
numerical  values  for  quality  characteristics  and  the  measure  of  limiting 
values  which  will  govern  for  purposes  of  the  criteria  should  be  defined. 
The  increasing  importance  of  this  concept  will  be  felt  as  and  when  the 
stream  quality  standards  are  implemented  and  enforced,  particularly  with 
the  advent  of  growing  use  of  automatic  stream  monitoring  installations. 


SCOPE  OF  THE  STUDY 

The  present  study  is  limited  to  the  following: 

1.  The  relative  importance  of  the  variations  in  reaction 
velocity  coefficients  on  the  dissolved  oxygen  response  in  the  stream 
are  evaluated  using  a  hypothetical  stream  situation  based  on  Streeter- 
Phelps*  formulation. 

2.  The  Ohio  River-Cincinnati  Pool  survey  data  for  the  river 
reach  between  miles  474.6  and  479.05,  published  by  the  U.  S.  Department 
of  Health,  Education  and  Welfare  (1960)  are  used  for  determining  the 
values  of  the  velocity  coefficients  affecting  deoxygenation  rates  in 
river  samples. 

3.  The  data  collected  by  the  Tennessee  Valley  Authority 
(1962)  for  the  prediction  of  stream  reaeration  rates  are  used  for 
formulating  the  hypothesis  concerning  the  chance  variation  of  the 
velocity  coefficient  governing  the  rate  of  atmospheric  reaeration 
in  rivers. 

4.  Statistical  models  describing  the  chance  variations  of 
the  two  velocity  coefficients  mentioned  above  are  formulated  and  tested. 

5.  Using  Monte  Carlo  techniques,  two  sets  of  data  for  the 
velocity  coefficients  are  generated  on  the  basis  of  the  statistical 
models  formulated.  The  generated  data  are  tested  and  verified. 

6.  The  generated  data  are  used  in  the  conceptual  model  for 
the  dissolved  oxygen-biological  oxygen  demand  (DO-BOD)  relationship  in 
streams  for  predicting  the  DO  response.  The  conceptual  model  is  verified 
using  the  Ohio  River-Cincinnati  Pool  survey  data  and  the  probabilistic 
variation  of  the  DO  responses  are  studied. 


II.   LITERATURE  REVIEW 

The  philosophy  of  protecting  the  stream  for  most  uses,  even 
future  uses,  is  being  increasingly  accepted,  since  it  is  good  from  the 
standpoint  of  the  water  user  and  the  public*  It  permits  planned  growth 
of  water-using  industries  and  recreation  side  by  side,  and  it  prevents 
undue  time  lags  between  the  need  for  water  use  and  its  realization.  In 
recent  years,  there  has  been  a  great  deal  of  attention  towards  the 
systems  analysis  approach  to  water  quality  management  considering  the 
river  flow  together  with  the  possibilities  of  affecting  or  controlling 
the  quality  of  the  water  at  various  use  points  by  dams,  water  purifi- 
cation plants,  and  waste  water  and  industrial  waste  treatment  plants 
(Thomas  and  Burden,  1963;  Liebman,  1965;  Loucks,  1965;  Montgomery, 
1964;  Worley,  1963).  In  most  of  these  studies  (Loucks,  1965;  Liebman, 
1965;  Montgomery,  1964;  Worley,  1963),  the  critical  dissolved  oxygen 
resulting  from  the  addition  of  organic  waste  to  a  water  course  is  taken 
as  a  controlling  parameter  for  studying  the  system  performance,  where 
critical  dissolved  oxygen  concentration  is  defined  as  the  minimum  con- 
centration of  dissolved  oxygen  in  the  river  below  a  waste  outfall.  If 
such  an  approach  is  to  yield  more  reliable  information,  it  is  imperative 
that  the  DO  response  of  the  receiving  stream  for  the  waste  loads  placed 
in  it  must  be  predicted  with  a  greater  degree  of  accuracy  than  is  pos- 
sible with  the  methods  adopted  in  these  studies. 

The  wastes  discharged  into  a  stream  from  municipal  and  in- 
dustrial treatment  plants  contain  a  large  variety  of  chemical  compounds. 
Of  primary  interest  in  this  work  is  that  portion  of  the  waste  which  is 
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biodegradable  and  henoe  oxygen  consuming*  When  this  material  is 
introduced  in  a  water  course,  it  undergoes  biochemical  oxidation, 
caused  by  microorganisms  which  utilise  the  organic  matter  for  energy 
and  growth. 

If  sufficient  dissolved  oxygen  is  present  in  the  water,  the 
microflora  are  primarily  aerobic,  uti listing  the  dissolved  oxygen  to 
carry  out  oxidation  reactions  producing  water  and  carbon  dioxide  as 
end  products.  If,  however,  sufficient  oxygen  is  not  present,  anaerobic 
organisms  predominate  resulting  in  undesirable  end  products.  Also  the 
concentration  levels  of  oxygen  have  profound  effects  on  the  physiology 
and  type  of  fishes  found  in  the  water  bodies.  The  lethal  effect  of 
low  concentrations  of  DO  appears  to  be  increased  by  the  presence  of 
toxic  substances  such  as  excessive  dissolved  carbon  dioxide,  ammonia, 
cyanides,  zinc,  lead  or  copper. 

In  situations  where  maintenance  of  well  balanced  fish  habitat 
is  one  of  the  primary  objectives,  the  stipulated  standards  for  dis- 
solved oxygen  are  that  it  shall  not  be  less  than  5.0  mg/1  during  at 
least  16  hours  of  any  24-hour  period,  nor  less  than  3,0  mg/1  at  any 
time  (Illinois  Sanitary  Water  Board,  1966;  Aquatic  Life  Advisory  Com- 
mittee of  ORSANCO,  1960),  Where  other  uses  like  the  industrial  water 
are  more  important,  the  stream  standards  adopted  for  such  uses  with 
respect  to  DO  differ  from  the  above. 

Because  of  the  variety  of  oxygen -demanding  organ! cs  in  the 
wastes,  it  is  common  to  measure  the  strength  of  wastes  in  terms  of 
their  biochemical  oxygen  demand  (BOD)  rather  than  to  analyze  for  the 
chemical  constituents  of  the  wastes. 
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The  most  widely  used  mathematical  models  for  predicting  the 
DO  responses  in  a  stream  are  either  the  one  proposed  by  Streeter  and 
Phelps  (1925)  or  the  modified  predictor  equations  based  on  Streeter- 
Phelps'  formulation  (Camp,  1963;  Dobbins,  1964),  though  a  few  other  con- 
cepts in  DO  predictions  (Churchill  and  Buckingham,  1956;  Thayer  and 
Krutchkoff,  1965)  have  been  put  forth  and  applied  with  varying  degrees 
of  success* 

The  Streeter -Phelps  equation  considers  only  two  mechanisms 
affecting  DO,  namely  the  removal  of  oxygen  by  bacterial  oxidation  of 
organic  matter  and  the  absorption  of  oxygen  from  the  atmosphere.  First 
order  kinetics  are  employed  to  express  deoxygenation  and  reoxygenation. 
The  rate  of  removal  of  BOD  on  the  basis  of  first  order  kinetics  is 
given  by: 


i  -  -  h1  (l) 


where  L  is  the  ultimate  first  stage  BOD  (milligram  per  liter)  remaining 
to  be  satisfied  at  any  time  t  (days)  and  K^  (base  e)  is  the  reaction 
rate  coefficient,  which  depends  on  the  characteristics  of  the  waste, 
temperature,  type  of  microorganisms  present,  and  other  environmental 
factors.  The  integration  of  this  equation  yields 

-Klt 
L  -  Lae  (2) 

where  La  is  the  initial  (t  «  0)  ultimate  first  stage  demand*  Since  BOD 

is  measured  in  terms  of  oxygen  consumed,  the  rate  of  oxygen  consumption 

equals  the  rate  at  which  BOD  is  satisfied,  i.e.,  $pL  «  -  ^  -  K*L,  where 

dt     dt    x 
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D  is  the  dissolved  oxygen  deficit  defined  as  the  difference  between 
saturation  concentration,  Cg,  at  the  river  water  temperature  and  the 
oxygen  concentration,  C,  in  the  river. 

Reoxygenation  process  in  the  stream  due  to  atmospheric  re- 
aeration  is  also  considered  as  a  first  order  reaction  which  depends  on 
the  dissolved  oxygen  deficit,  D,  This  can  be  expressed  mathematically 
as: 


£  -  K2  (C,  -  C)  -  K2D  -  -  §  (3) 


where  K2  (base  e)  is  the  reaction  rate  coefficient  defining  the  re- 
aeration  process* 

Combining  the  rates  of  these  two  reactions,  and  writing  the 
resulting  equation  in  terms  of  dissolved  oxygen  deficit,  the  response 
of  the  receiving  stream  for  a  single  waste  source  is  defined  by  the 
differential  equation  as: 


f  -  KlL  .  K2D  (4) 


The  solution  of  this  differential  equation  with  appropriate  initial 
conditions  yields  the  well  known  and  classical  Streeter-Phelps  equation 
for  stream  self  purification  capacity  and  is  given  as: 

KlLa 
D  "  k2  -  Kx  0*P  <-Klfc)  -  «P  (-K2O]  +  Da  e*P  <-K2t)   (5) 

where  Da  is  the  initial  dissolved  oxygen  deficit  in  the  river  reach  at 
time  t  «  0  and  all  other  terms  are  as  defined  previously. 
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Figure  1  shows  the  typical  oxygen  sag  curve.  The  point  of 
maximum  deficit  (or  minimum  dissolved  oxygen  concentration)  is  known 
as  the  critical  point*  The  critical  time,  or  time  at  which  the  maximum 
deficit  occurs  (tc),  and  the  corresponding  maximum  or  critical  deficit 
(Dc)  are  defined  by  the  following  equations  as: 

*        K2      (K2  -  ^i' 
<=c  -  K7TI7  {m  Ck^  (1  -   KlLa*  D.)]}        (6) 

and 

Kl 
Dc  "  kJ  La  exP  <"*!*<:>  <7> 

In  the  river  reach  prior  to  the  critical  point,  the  rate  of  deoxygena- 
tion  is  greater  than  the  reaeration  rate  and  beyond  this  point,  the 
converse  holds.  The  river  is  said  to  recover  in  the  latter  region  from 
the  pollutional  load  discharged  into  it* 

The  shortcomings  and  criticisms  expressed  against  the  Streeter- 
Phelps  equation,  Eq*  5,  for  stream  self  purification  capacities  and 
attempts  of  other  investigators  to  modify  these  concepts  in  order  to 
improve  upon  the  predictions  are  discussed  in  the  following  paragraphs. 
There  are  also  altogether  different  approaches  taken  for  the  prediction 
of  DO  response  of  the  streams  to  waste  loads  placed  in  it  and  these  are 
discussed  along  with  the  suggested  modifications  for  the  Streeter-Phelps 
concepts  c 

The  Streeter-Phelps  sag  equation  when  fitted  to  stream  DO  and 
BOD  data  have  yielded  such  a  wide  range  of  values  of  the  parameters  Kj 
and  K«  as  to  suggest  that  some  important  factors  have  not  been  taken 


12 


•—I  c      -— 

^^ 

B 

mm* 

u 
m 
u 
u 

s 


o 

(8 
CI 

25 


T 

Da 

* 

; 

, 

D                ^"^ 

— 

D 

c 

*'yr 

« tc 

* 

* 

—  t   

H 

_______ 

Tine,   days 


FIGURE   1.      DISSOLVED  OXYGEN  SAG  CURVE 
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into  account  in  the  formulation  or  that  the  variables  have  not  been 
incorporated  in  the  correct  mathematical  form.  Several  authors  (Camp, 
1963;  Dobbins,  1964;  O'Connor,  1967)  have  proposed  modifications  to 
the  Streeter-P helps  formulation  by  taking  into  account  other  parameters 
which  affect  the  oxygen  uptake  in  the  stream  like  removal  of  BOD  due 
to  adsorption  and  sedimentation,  addition  of  BOD  from  benthal  layers 
to  the  overlying  water,  benthal  oxygen  demand,  oxygen  addition  due  to 
photosynthesis,  algal  respiration,  etc. 

Some  apprehension  also  exists  concerning  the  validity  of  the 
application  of  the  first  order  concept  and  how  the  results  of  the  BOD 
test  should  be  interpreted  mathematically  (Gates,  1966;  Thomas,  1961; 
Young,  1965).  Thomas  (1961)  assumes  second  order  (bimolecular)  reaction 
kinetics  for  BOD  removal,  in  his  "Step  Method"  for  estimating  the  oxygen 
uptake  in  a  stream.  The  advantage  claimed  by  the  author  is  that  the 
variation  of  the  BOD  reaction  velocity  parameter  is  considerably  re- 
duced if  a  second  order  rate  equation  is  used  in  piece  of  the  first 
order  equation.  Young  (1965)  concludes,  based  on  the  analysis  of  the 
data  for  sewage,  that  the  first  and  second  order  kinetics  fit  the  first 
stage  BOD  data  with  equal  precision  at  20°C  and  35°C.  Since  the  validity 
of  the  second  order  kinetics  concept  over  the  first  order  kinetics  has 
not  been  well  established,  most  of  the  recent  works  (Liebman,  1965; 
Loucks,  1965;  O'Connor,  1967;  Thayer,  1966)  use  first  order  kinetics 
for  the  oxygen  depletion  due  to  bacterial  degradation  of  organic 
matter.  Inasmuch  as  this  concept  is  more  widely  recognized  and  used 
in  the  sanitary  engineering  field,  the  BOD  progression  will  be  assumed 
to  follow  first  order  kinetics. 
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The  merit  of  employing  a  mathematical  model  developed  from 

data  obtained  in  such  a  nondynamic  system  as  BOD  bottle  and  applying 

it  directly  to  a  dynamic  system  such  as  a  stream  has  been  questioned 

in  the  recent  past  (Gates,  1966;  Gannon,  1963;  Nejedly,  1966),  Gannon 

(1963)  states 

•••it  is  apparent  that  no  adequate  explanation  is 
available  to  account  for  the  differences  between 
the  laboratory  kj  and  the  river  kr.  This,  no 
doubt,  is  partially  due  to  inadequacies  of  the 
existing  BOD  techniques.  Probably  it  may  never 
be  possible  to  develop  a  BOD  procedure  in  the 
laboratory  that  will  duplicate  the  natural  river 
conditions,  but  certainly  efforts  should  be 
directed  toward  improving  existing  procedures. 

In  addition  to  improvement  of  the  laboratory 
BOD  test  it  appears  that  much  useful  fundamental 
information  could  be  obtained  about  BOD  removal 
rates  from  controlled  studies  using  artificial 
channels.  This  approach,  then,  would  be  somewhere 
between  the  natural  stream  conditions  and  the 
highly  artificial  bottle  methods, 

Isaacs  and  Gaudy  (1967)  from  their  studies  specifically 
designed  to  gain  an  insight  into  the  equivalence  of  K^  determined 
in  the  bottle  and  in  a  laboratory  stream  model,  using  synthetic 
sewage  consisting  of  equal  parts  of  glucose  and  glutamic  acid  came 
to  the  conclusion  that  data  obtained  in  BOD  bottles  could  be  used  to 
make  a  reasonable  prediction  of  the  critical  DO  provided  the  seed 
concentration  in  the  BOD  bottle  was  identical,  or  nearly  identical, 
to  the  biological  solids  concentration  in  the  receiving  stream.  They 
further  observe,  based  on  the  model  stream  study,  that  a  fairly  ac- 
curate prediction  could  be  made  by  estimating  the  extent  of  the  first 
stage  exertion*  Since  these  conditions  could  be  easily  met  by  making 
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observations  on  river  samples  obtained  downstream  of  a  waste  source, 
the  idea  of  extrapolating  the  BOD  bottle  data  to  river  situations 
cannot  be  ruled  out,  at  least  until  a  better  method  could  be  devised 
and  established* 

Nejedly  (1966)  holds  the  view  that  longitudinal  mixing,  i.e., 
mixing  in  the  direction  of  flow,  induces  contacts  between  particles  of 
decaying  organic  substances  with  different  detention  periods  within  the 
river  stretch  considered,  and  hence  a  more  rapid  BOD  reduction  takes 
place  in  a  stream  than  in  a  bottle.  This  view  is  opposed  to  those  of 
Camp  (1965),  Dobbins  (1964),  0,Connor  (1967)  and  several  others  who 
hold  that  the  effect  of  longitudinal  dispersion  is  negligible  in  deter- 
mining streams*  assimilative  capacity*  Though,  Nejedly  presents  data 
to  support  his  hypothesis,  no  quantitative  conclusion  for  the  effect  of 
longitudinal  mixing  on  the  reaction  coefficient  K  could  be  determined. 

It  is  generally  recognized  that  the  differences  in  removal 
rates  of  BOD  in  BOD  bottles  and  in  actual  streams  are  due  to  the  phen- 
omena of  biological  accumulation  and  removal  of  BOD  due  to  sedimentation 
without  the  DO  being  affected  (Gannon,  1963;  O'Connor,  1967).  These  two 
removal  rates  are  taken  as  the  same  in  the  absence  of  settleable  or- 
ganic solids  and  fixed  aquatic  vegetations  in  the  receiving  stream 
(O'Connor,  1967). 

In  spite  of  the  fact  that  the  bottle  BOD  technique  has  been 
subjected  to  severe  criticism,  it  has  been  used  in  one  form  or  the 
other  in  determining  stream  assimilative  capacities.  Recognizing  that 
deoxygenation  due  to  bacterial  degradation  of  organic  matter  is  only 
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one  of  the  several  mechanisms  which  affect  the  BOD  removal  in  a 
stream,  there  seems  to  be  ample  justification  in  asing  the  data 
obtained  in  BOD  bottles  to  determine  river  self  purification  capac- 
ities, particularly  when  the  mechanisms  affecting  the  BOD  removal 
without  oxygen  consumption  are  considered  separately  in  mathematical 
modeling* 

Another  severe  criticism  voiced  against  the  use  of  Streeter- 
Phelps  equation  or  its  modifications  is  that  the  so-called  reaction 
velocity  constants  Kj,  K£,  etc,  are  found  to  vary  considerably. 
LeBosquet  and  Tsivoglou  (1950)  observe: 

The  velocity  constant  for  reoxygenation,  K£,  is  subject 
to  wide  variation  in  different  streams,  in  individual 
stretches  of  any  one  stream  and  for  various  river 
stages  in  the  same  stretch.  Although  variations  in 
K^  are  usually  of  somewhat  lesser  magnitude,  they 
too  may  be  appreciable,  and  in  many  situations  it 
is  quite  difficult,  if  not  impossible,  to  derive 
acceptable  values  for  these  velocity  coefficients. 

Churchill  and  Buckingham  (1956)  state  regarding  the  variability  of 

these  coefficients  as  follows: 

The  so  called  constants  K^  and  K2  are  determined 
from  an  intensive  survey,  usually  at  a  fairly  low, 
steady,  flow.  However,  if  another  survey  is  made  on 
the  same  reach,  even  fairly  soon  after  the  first 
one,  it  at  once  becomes  obvious  that  the  "constants** 
may  vary  considerably  from  the  values  first  deter- 
mined ....  Therefore,  when  the  derived  "constants" 
are  far  from  constant,  the  results  of  computation  in- 
volving these  "constants"  for  other  conditions  of 
load,  flow,  and  temperature,  are  certainly  open  to 
question. 

The  main  aim  of  this  work  as  stated  earlier  is  to  ascertain  the  sig- 
nificance of  the  variations  in  the  velocity  coefficients  K^  and  K~  in 
defining  the  DO  response  of  a  receiving  stream  and  to  develop  a  procedure 


17 

for  determining  the  DO  taking  the  variability  in  these  rate  coefficients 
into  consideration,  if  these  variations  are  significant. 

Velz  (1939)  presented  a  method  for  oxygen  balance  studies  in 
polluted  rivers  in  which  the  atmospheric  river  reaeration  is  computed 
from  an  empirical  relationship,  established  in  a  quiescent  column  of 
water,  defining  the  relationship  of  the  terms:  percent  saturation 
after  a  known  elapsed  time,  linear  depth  of  column,  molecular  diffusion 
coefficient,  and  the  initial  percent  saturation  of  the  column  of  water. 
The  application  of  an  empirical  relationship  established  in  a  quiescent 
body  of  water  to  a  dynamic  system  as  in  a  river  situation  is  highly 
questionable,  particularly  when  the  effect  of  molecular  diffusion  of 
oxygen  is  considered  negligible  in  a  stream.  Also  a  family  of  curves, 
giving  the  relationship  between  BOD  and  time  at  temperatures  0°C  to 
30°C  are  presented  for  the  purpose  of  determining  the  BOD  reduction 
along  the  course  of  the  river.  Though  first  order  reaction  kinetics  is 
assumed  in  defining  these  curves,  no  distinction  is  made  for  the  differ- 
ences in  the  characteristics  of  the  wastes  under  examination,  i.e.,  the 
percent  reduction  of  the  pollutional  load  in  a  given  flow-through  time, 
is  taken  as  independent  of  the  nature  of  the  waste.  Velz's  method 
suffers  from  the  same  shortcomings  discussed  earlier  for  the  Streeter- 
Phelps  formulation  and  in  addition  this  method  does  not  take  cognizance 
of  the  nature  of  wastes  which  is  an  important  factor  in  biodegradability. 

Thomas  (1948)  discussed  a  technique  for  stream  analysis  which 
recognizes  that  the  K,  value  may  not  be  0.23,  the  commonly  accepted 
value  in  sanitary  engineering  practice,  and  also  that  there  may  be  BOD 
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sinks  other  than. that  represented  by  bio-oxidation.  Furthermore,  his 
discussion  indicates  that  he  did  not  consider  these  additional  BOO 
sinks  to  be  oxygen  sinks.  This  technique  requires  that,  through  the 
use  of  the  BOD  test,  the  values  of  K^  and  ultimate  first  stage  demands 
be  determined  at  boundary  stations  of  a  hydraulically  uniform  stretch 
of  stream.  The  difference  between  the  ultimate  first  stage  BOD  values 
of  the  upstream  and  downstream  stations  represents  the  reduction  in 
oxygen  demand  realized  through  the  reach.  This  information  combined 
with  the  time  of  passage,  and  assuming  first  order  kinetics,  is  used 
to  determine  the  rate  constant  expressing  the  reduction  in  potential 
oxygen  demand  by  the  biomass  is  determined  by  averaging  the  K^  values 
observed  for  samples  from  the  two  boundary  stations.  The  difference 
between  this  rate  and  the  rate  observed  for  the  stream  represents  the 
rate  at  which  potential  oxygen  demanding  material  is  going  to  non- 
oxygen  consuming  sinks. 

Although  this  approach  recognizes  that  Kj  may  vary,  it  assumes 
that  the  stream  values  of  K.  and  L  can  be  predicted  from  bottle  measure- 
ments. The  concept  that  not  all  potential  oxygen  consumption  need 
necessarily  be  satisfied  by  actual  consumption  was  a  significant  con- 
tribution to  stream  analysis.  But  for  this  modification,  Thomas'  tech- 
nique is  identical  to  that  employed  by  Streeter  and  Phelps  (1925)  and 
subject  to  the  same  previously-indicated  reservations. 

LeBosquet  and  Tsivoglou  (1950)  published  data  for  the  Ohio 
River  indicating  a  radically  different  method  of  pollutional  analysis 
employing  statistical  methods.  They  established  a  linear  relationship 
between  critical  DO  deficit  in  mg/1  and  reciprocal  of  river  flow  in 
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cubic  feet  per  second,  with  a  high  degree  of  correlation.  This  method 
permits  direct  analysis  of  oxygen-flow  relationship  in  a  given  stretch 
avoiding  the  use  of  factors  such  as  Kj  and  K2.  The  method  outlined  by 
them  cannot  give  better  results  than  the  oxygen  sag  formula  since  the 
characteristics  of  the  wastes  and  the  reoxygenation  capacities  of  the 
stream  fluctuate  greatly  and  such  variations  will  introduce  inaccuracies. 
It  is  based  mainly  on  the  surmise  that  the  effect  of  the  variations  in 
K^  and  K2  are  relatively  minor,  which  is  not  a  valid  conjecture  as  will 
be  borne  out  by  this  study. 

Churchill  and  Buckingham  (1956)  extended  the  concept  of 
LeBosquet  and  Tsivoglou  (1950)  and  presented  an  analysis  employing 
statistical  techniques  of  multiple  correlation  and  giving  the  maximum 
deficit  in  terms  of  a  linear  equation  with  5-day  BOD,  water  temperature 
and  stream  discharge  as  variables.  The  advantages  claimed  for  this 
method  are  that  evaluation  of  factors  like  K^,  K2,  L  and  flow-through 
time  are  avoided  and  that  the  correlation  analysis  predicts  the  critical 
value  reasonably  well. 

The  authors  (Churchill  and  Buckingham,  1956)  recommend 
choosing  a  station  for  making  5-day  BOD  observations  downstream  of 
the  pollution,  and  as  close  to  the  critical  section  as  possible.  Since 
the  critical  section  in  the  river  fluctuates  with  variations  in  stream 
flow,  and  the  characteristics  of  the  waste,  the  criteria  suggested  by 
the  authors  for  selection  of  site  for  making  BOD  observations  to  be 
used  in  the  regression  analysis  is  not  easily  met  since  it  is  not 
feasible  to  estimate  the  location  of  the  critical  section  from  their 
statistical  model.  Also,  it  is  not  possible  to  predict  the  state  of 
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oxygen  concentration  in  different  sections  of  the  river  using  the 
correlation  technique*  It  is  of  paramount  importance  to  estimate 
assimilative  capacity  at  any  desired  reach  and  not  merely  the  critical 
reach  of  the  river,  particularly  when  comprehensive  river  basins 
planning  are  envisaged  with  planned  development  of  rivers. 

Thomas  (1961)  presented  a  method  called  the  "Step  Method"  for 
estimating  oxygen  uptake  in  a  stream.  Like  the  St  reeter-P  helps  formu- 
lation, only  two  mechanisms  affecting  DO  in  a  stream  are  considered, 
namely  the  removal  of  oxygen  due  to  bacterial  decomposition  and  addition 
of  oxygen  from  the  atmosphere,  even  though  he  has  emphasised  earlier 
(1948)  the  importance  of  considering  removal  and  addition  of  BOD  due  to 
sedimentation  and  scour  respectively.  The  only  departures  from  the 
classical  Streeter-Phelps  formulation  are  that  the  bacterial  degradation 
of  organic  matter  is  assumed  to  be  governed  by  second  order  kinetics 
and  that  the  values  of  5-day  BOD  are  used  instead  of  the  ultimate  first 
stage  BOD.  In  the  light  of  the  investigations  of  Young  and  Clark 
(1965).  and  Isaacs  and  Gaudy  (1967)  discussed  earlier,  the  "Step  Method" 
does  not  appear  to  improve  the  predictive  techniques  for  stream  self 
purification  capacities. 

Camp  (1963).  and  Dobbins  (1964)  proposed  equations  which  are 
very  similar,  for  predicting  dissolved  oxygen  deficit  by  taking  into 
account  benthal  demand,  removal  or  addition  of  BOD  due  to  deposition  or 
scour,  and  photosynthesis.  0,Connor  (1967)  further  extended  these  ideas 
by  considering  the  diurnal  variation  of  oxygen  addition  due  to  photosyn- 
thesis instead  of  treating  it  as  time  constant  for  the  whole  24-hour 
period  as  in  Camp's  model,  and  adding  another  term  for  algal  respiration 
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in  the  oxygen  sinks.  Since  these  formulations  (Camp,  1963;  Dobbins, 
1964;  O'Connor,  1967)  are  only  modifications  of  the  basic  Streeter- 
Phelps  formulation,  considering  a  few  more  plausible  mechanisms 
affecting  the  DO-BOD  relationship  in  streams,  they  do  not  obviate 
all  the  criticisms  put  forth  for  the  Street«r-P helps  equation. 

Thomann  (1963)  postulated  a  mathematical  model  using  the 
systems  analysis  concept  for  predicting  the  response  of  any  water 
quality  parameter,  affected  either  by  conservative  or  nonconservative 
pollutants.  The  model  is  developed  for  a  one-dimensional  estuary  for 
which  stream  analysis  is  only  a  particular  case.  The  model  assumes 
that  a  body  of  water  can  be  segmented  into  a  discrete  number  of  parts 
each  of  which  is  treated  as  a  stationary  linear  subsystem  forming  part 
of  the  whole  system.  The  time  rate  of  change  of  a  water  quality  parameter 
particularly  DO,  in  any  of  the  subsystems  is  taken  as  affected  by  the 
advection  and  diffusion  of  DO  both  into  and  out  of  the  system  and  by 
the  sources  and  sinks  of  DO  within  each  subsystem.  These  additions  or 
removals  of  DO  are  considered  as  the  variables  that  force  the  DO  to 
respond.  An  equation  representing  statinary  linear  subsystem  for  DO 
concentration  at  a  given  time,  with  the  sources  and  sinks  of  DO  as 
forcing  functions  and  also  with  feed-back  mechanisms  from  adjoining 
subsystems,  is  developed.  Since  the  subsystems  that  form  the  overall 
system  are  all  linear,  the  solutions  obtained  for  each  subsystem  are 
linearly  superimposed  to  predict  the  response  of  the  system  as  a  whole. 

Though  Thomann  (1963)  has  indicated  conceptually  a  method  for 
predicting  the  DO  response  when  the  decay  coefficient  for  organic  matter 
changes  with  time,  no  mention  is  made  in  his  theoretical  development, 
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how  the  variations  with  tiiae  in  the  decay  coefficients  are  evaluated* 
Moreover,  the  possible  variations  in  the  reaeration  coefficient  are 
overlooked.  It  is  well  documented  (TVA,  1962)  that  the  maximum  re- 
aeration  coefficient  observed  under  constant  flow  conditions  in  the 
Tennessee  Valley  rivers,  when  the  temperatures  remained  nearly  uniform, 
ranged  approximately  from  2  to  9  times  the  minimum  values  depending  on 
the  characteristics  of  the  river,  Thomann's  work  marks  a  significant 
advance  in  pollution  analysis  wherein  an  attempt  has  been  made  to  con- 
sider all  the  factors  affecting  DO  in  a  body  of  water  taking  into  account 
the  possible  variations  in  some  of  the  parameters  affecting  the  system. 

Thayer  (1966)  viewed  the  oxygen  relationship  in  streams  as  a 
stochastic  birth  and  death  process  with  the  BOD  and  DO  being  increased 
and  decreased  by  small  amounts  in  a  very  short  interval  of  time.  The 
stochastic  model  provides  for  the  joint  density  function  for  both  pol- 
lution and  dissolved  oxygen  for  different  initial  conditions.  The  ad- 
vantage claimed  by  the  author  is  that  in  addition  to  predicting  the  mean 
DO  concentration  which  is  the  same  as  that  predicted  by  the  determin- 
istic equation  of  Dobbins  (1964),  it  affords  a  measure  of  the  variance 
of  DO  from  its  mean  value,  thus  enabling  one  to  predict  the  probability 
of  DO  downstream  of  a  source  of  pollution  taking  any  given  value. 

In  developing  his  model,  Thayer  (1966)  did  not  consider  the 
effects  of  photosynthesis  and  algal  respiration.  The  prediction  equation 
was  verified  with  the  data  published  by  the  Resources  Agency  of  California 
(1962)  for  the  Sacramento  River.  Though  the  published  data  indicate  a 
value  for  K.  ranging  from  0.12  to  0.55  and  emphasise  the  importance  of 
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photosynthesis  and  algal  respiration  in  the  river  reach  below  Sacra- 
mento, these  effects  were  totally  ignored  in  the  model  verification* 
Also  in  the  controlled  laboratory  experiments  with  dextrose  as  the 
only  substrate  for  simulating  river  conditions,  values  ranging  from 
0.076  to  0.167  were  obtained  forKj  (Thayer,  1966).  Whereas,  in  the 
verification  of  his  theory  with  the  results  of  the  laboratory  data,  a 
value  of  0.07  for  Kj  has  been  assumed  on  the  supposition  that  the  as- 
sumed value  is  reasonably  close  to  the  observed  values.   In  the  opinion 
of  this  writer,  the  verification  of  the  mathematical  model  (Thayer, 
1966)  either  with  the  published  data  for  the  Sacramento  River  or  with 
the  laboratory  data  leave  much  to  be  desired. 

From  the  foregoing  general  discussion  on  the  present  state 
of  knowledge  concerning  the  stream  assimilative  capacities,  it  can  be 
summarised  that  though  there  are  strong  indications  from  field  measure- 
ments for  the  variability  in  the  parameters  affecting  DO-BOD  relation- 
ship, and  though  a  need  for  talcing  into  account  these  variations  is 
strongly  felt,  all  the  theories  so  far  postulated  except  that  of  Thomann 
(1963)  fail  to  consider  this  aspect.  Even  though  there  is  considerable 
criticism  for  extrapolating  the  BOD  bottle  observations  to  river  con- 
ditions, this  is  the  only  expedient  way  available  at  the  present  time 
for  defining  the  biological  interactions  in  the  environment.  The  con- 
cept of  first  order  kinetics,  first  introduced  by  Streeter  and  Phelps 
(1925),  with  the  parameters  determined  by  bottle  BOD  tests  will  continue 
to  be  used  until  a  better  and  more  reliable  method  is  evolved  to  sub- 
stitute this  concept. 
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III.   SENSITIVITY  ANALYSIS 

DEFINITION 

Sensitivity  analysis  applies  to  the  concept  in  which  all 
except  one  of  the  variables  in  a  system  are  held  constant  and  this 
single  exception  is  allowed  to  vary  through  its  full  possible  range. 
The  effects  on  certain  system  performance  are  noted,  testing  the 
sensitivity  of  the  system  to  that  particular  variable  which  is  varied. 
If  the  measured  output  varies  significantly,  the  system  is  said  to 
be  sensitive  to  that  variable  for  the  given  conditions.   In  a  mathe- 
matical sense,  the  slope  of  the  output  as  a  function  of  input  is  a 
sensitivity  measure. 

Two  variables  of  a  system  may  be  permitted  to  vary  simultan- 
eously, each  having  an  effect  on  the  output  measure.  Mathematically 
this  becomes  a  three-dimensional  model;  the  output  function  becomes  a 
surface  rather  than  a  line.  This  is  called  a  response  surface.  A  well 
defined  response  surface  contains  outcomes  of  all  possible  combinations 
of  the  two  variables.   In  a  mathematical  sense,  the  sensitivity  of 
the  response  surface  to  a  single  independent  variable  is  a  partial 
derivative  of  the  response  surface  equation  with  respect  to  that  variable. 

PROCEDURE 

The  suggested  general  procedure  for  carrying  out  a  sensitivity 
analysis  is  as  follows: 

1.  The  independent  variables  of  the  system  are  listed  and 
the  possible  range  of  values  for  each  of  the  variables  is  ascertained. 
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2.  The  value  of  the  response  surface  variable  is  calculated 
using  nominal  values.  This  value  serves  as  a  reference  point* 

3.  Percent  change  in  the  response  surface  is  computed  by 
substituting  the  low  value  of  one  of  the  independent  variables, 

4.  The  low  value  is  replaced  with  an  incremented  value  for 
the  same  variable  as  in  step  3  and  the  percent  change  in  the  response 
surface  is  evaluated.  This  is  repeated  until  all  the  values  within 
the  possible  range  of  values  for  the  variable  under  consideration  are 
covered. 

5.  Steps  3  and  4  are  repeated  for  each  of  the  independent 
variables  which  is  likely  to  assume  a  range  of  values  instead  of  a 
single  value. 

The  above  procedure  for  individual  variables  analysis  can  be  modified 
to  test  for  changes  of  variables  in  groups  which  may  be  related. 

Nominal  Values  in  the  Streeter -Phelps  Equation 

The  Streeter -Phelps  equation  for  stream  assimilative  capacity 
(Eq.  5)  will  be  used  in  the  sensttivity  analysis  for  determining  the 
significance  of  the  variations  in  reaction  velocity  coefficients  K.  and 
Kg  in  predicting  the  dissolved  oxygen  response.  It  is  generally  ac- 
cepted in  sanitary  engineering  practice  that  the  average  value  for  K^ 
at  20°C  for  domestic  wastes  is  0.23  per  day  (base  e)  (Fair  and  Geyer, 
1954;  Sawyer,  1960),  even  though  values  ranging  from  0.1  to  0,7  per  day 
are  indicated  by  some  authors  for  domestic  and  industrial  wastes  (Camp, 
1963;  Eckenf elder  and  O'Connor,  1961).  The  Kj^  value  at  20°C  reported 
for  river  samples  collected  downstream  of  Sacramento  sewage  treatment 
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plant  discharge  is  found  to  range  from  0,12  to  0.55  per  day  with  an 
average  value  of  0.294  per  day  (The  Resources  Agency  of  California, 
1962).  The  computed  Kj  values  for  the  Ohio  River  samples  obtained 
during  a  1957  survey  (USPHS,  1960)  in  the  river  reach  between  miles 
474.6  and  479,05,  where  all  the  pollution  entered  prior  to  the  reach 
under  consideration,  showed  a  range  from  0.05  to  0.32  per  day  with  an 
average  value  of  0.173  per  day.  From  the  average  value  cited  for 
K.  in  the  literature  and  from  the  values  of  K  obtained  in  actual 
river  surveys,  it  can  be  concluded  that  the  nominal  (average)  value 
for  the  purpose  of  sensitivity  analysis  can  be  taken  as  0.23  per  day 
at  20°C. 

The  nominal  value  for  K2  at  20°C  is  not  very  well  defined 
since  it  depends  on  the  characteristics  of  the  stream*  The  classifi- 
cation of  streams  from  reaeration  point  of  view  is  only  subjective  in 
nature.  Fair  and  Geyer  (1954)  indicate  for  large  streams  a  value  of  K2 
(base  e)  at  20°C  ranging  from  0,4  to  0.7  per  day.  Eckenfelder  and 
O'Connor  indicate  a  common  range  of  K2  from  0.2  to  10.0  per  day,  the 
lower  value  representing  deep  slow-moving  rivers  and  the  higher  value 
for  rapid  shallow  streams.   In  one  set  of  experiments  (No.  14)  in  the 
Holston  River,  which  could  be  classified  as  a  large  river  with  moderate 
velocity,  K2  had  a  range  of  values  from  0,10  to  1.18  per  day  with  an 
arithmetic  average  of  0.63  per  day  under  comparable  conditions  within 
practical  limits  (Buckingham,  1966),  Similar  ranges  of  values  were  ob- 
tained in  a  few  other  experiments  conducted  in  TVA  rivers.  Hence  an 
average  value,  in  the  case  of  large  streams,  for  K2  equal  to  0.6  per 
day  at  20°C  will  be  a  realistic  assumption. 
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Temperature  Effects  on  Reaction  Coefficients 

The  Importance  of  temperature  effects  on  the  velocity  coef- 
ficients K,  and  K2  has  long  been  veil  recognized  and  much  attention 
has  been  focused  on  this  aspect.  Streeter  and  Phelps  (1925)  described 
the  relationship  between  K^  and  temperature  by  the  expression 

K1(T)  "  K1<20)9T"20  (8> 

in  which  Kj^T)  is  the  deoxygenation  coefficient  at  any  temperature  T 
in  degrees  centigrade,  Ki(20)  *s  tne  deoxygenation  coefficient  at  20°C 
and  6  is  the  temperature  coefficient*  9  was  found  to  have  an  average 
value  of  1.047,  Theriault,  as  reported  by  Camp  (1963),  and  Thomas  (1961), 
confirmed  this  value  of  9,  based  on  his  Ohio  River  studies.  Recently, 
Zanoni  (1967)  in  an  attempt  to  establish  the  temperature  effects  on 
the  rate  of  deoxygenation  of  a  conventional  activated  sludge  waste 
water  treatment  plant  effluent,  found  the  value  of  9  to  be  1.048  for  the 
temperature  range  of  10  to  30°C.  Though  the  observed  values  for  9 
varied  within  a  small  range,  the  effects  of  this  variability  in  estimating 
the  deoxygenation  rate  at  other  temperatures,  knowing  its  value  at 
20°C,  will  not  be  significant.  This  temperature  relationship  with  9 
having  a  value  of  1.047  will  be  used  in  computing  the  nominal  values 
at  other  temperatures,  as  has  been  employed  in  earlier  works  (Ecken- 
f elder  and  O'Connor,  1961;  0*Connor,  1960;  Streeter  and  Phelps,  1925; 
Thomas,  1961;  Worley,  1963). 

Also,  it  has  long  been  known  from  experimental  evidence  that 
deoxygenated  water  will  absorb  oxygen  from  the  atmosphere  at  a  higher 
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rate  if  the  temperature  of  the  water  is  raised,  other  conditions  being 
held  constant.  In  the  normal  range  of  stream  temperatures,  a  rise  In 
water  temperature  results  in  a  decrease  in  viscosity,  density,  and 
surface  tension.  It  is  difficult  to  distinguish  the  exact  role  played 
by  each  of  these  factors.  Because  it  is  only  the  net  overall  effect 
of  temperature  on  re aeration  that  is  of  interest  in  most  investigations, 
experiments  have  been  conducted  in  the  past  for  measuring  reaeration 
rates  at  several  temperature  levels.  Results  of  such  experiments  have 
been  generally  reported  in  the  following  form: 

K2(T)  -  K2(20)*T"20  <*) 

in  which  K2(T)  is  the  reaeration  coefficient  at  temperature  T  in  degrees 
centigrade,  K2(20)  *s  the  reaeration  coefficient  at  20°C  and  $  is  the 
temperature  coefficient.  Several  authors  indicated  values  for  $  ranging 
from  1.008  to  1.02  as  reported  by  the  Committee  on  Sanitary  Engineering 
Research  in  their  Thirty-first  Progress  Report  (1961).  The  committee 
(1961)  concluded,  after  conducting  a  series  of  carefully  controlled 
experiments  that  the  thermal  coefficient  $  remains  constant  over  a 
wide  range  of  turbulence  conditions,  with  a  value  of  1.0241.  The  re- 
lationship shown  in  Eq.  9  with  $  equal  to  1.0241  will  be  used  in  com- 
puting the  nominal  values  at  other  temperatures. 

Table  1  shows  the  nominal  values  of  the  reaction  velocity 
coefficients  Kj  and  K2  at  temperatures  10°C,  20°C,  and  30°C  used  in 
the  sensitivity  analysis. 
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TABLE  1 
NOMINAL  VALUES  OF  REACTION  COEFFICIENTS  Kt  AND  K2 


Parameters 

Nominal  Values 

At  10°C 

At  20°C 

At  30°C 

Kl 
K2 

0.15 
0.47 

0.23 
0.60 

0.36 
0.76 

The  range  of  values  considered  for  Kj^  and  K2  are  0.1  to  0.6 
per  day  and  0.2  to  3,0  per  day  respectively  with  incremented  values 
of  0.1  per  day  for  Kj  and  0.2  per  day  for  K2.  The  initial  conditions 
assumed  for  ultimate  BOD  and  DO  deficit  in  the  hypothetical  stream  at 
the  source  of  pollution  are  20  mg/1  and  1.5  mg/1  respectively,  which 
are  considered  as  realistic  values  for  these  parameters,  having  been 
based  on  actual  observations  made  on  river  samples  (Gannon,  1963; 
USPH5,  1960).   Since  the  effect  of  temperature  on  ultimate  BOD  is 
considered  insignificant  (Gotaas,  1948;  Thomas,  1961;  Zanoni,  1967), 
no  correction  for  temperature  effects  will  be  made  for  estimating 
initial  ultimate  BOD  at  different  temperatures  considered  in  the 
sensitivity  analysis. 

RESULTS 

The  critical  DO  predicted  by  the  Streeter-Phelps  equation  is 
considered  for  evaluating  the  response  of  the  hypothetical  stream 
system.  The  value  obtained  for  critical  DO,  when  K^  and  K2  take  their 
respective  nominal  values,  is  treated  as  the  nominal  value  of  critical 
DO  for  the  given  conditions  of  initial  BOD,  DO  deficit,  and  temperature. 
Percent  deviations  of  the  system  responses  are  computed  when  K^  and  K2 
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take  different  values  within  the  possible  range  and  these  are  based  on 
the  reference  value  for  the  critical  DO  in  the  stream, 

Eq,  5,  which  expresses  the  relationship  between  D,  La,  Da»  t, 
K^  and  K2,  becomes  indeterminate  when  K^  and  K2  are  equal.   Since  the 
ranges  of  values,  which  these  two  coefficients  can  assume,  overlap  it 
is  necessary  to  predict  the  critical  deficit  when  Kj  equals  K2.  The 
solution  for  the  Streeter-Phelps  formulation,  Eq,  4,  is  given  below  for 
the  case  when  K^  and  K2  are  equal. 

D  -  KxLat  exp  (-l^t)  +  Da  exp  (-KLt)  (10) 

The  time  of  flow  for  critical  deficit,  and  critical  deficit 
are: 

La  -  Da 
'c  -  "kIET  (ID 

Dc  -  La  exp  (-Kxtc)  (12) 

A  computer  program  incorporating  the  general  Streeter-Phelps  equation, 
and  the  solution  for  the  situation  when  K^  equals  K2  was  written  in 
FASTRAN  langaage  for  the  IBM  7094-1401  system  to  evaluate  the  state  of 
oxygen  concentration  at  1,  2,  3,  4,  and  5  days  of  flow -through  time  in 
order  to  trace  the  oxygen  sag  profile.  The  program  computes  the  flow- 
through  time  for  critical  deficit  and  hence  the  critical  deficit  and 
critical  DO  for  different  combinations  of  values  for  K}  and  K2  with  the 
assumed  initial  conditions  for  La  and  Da.  The  flow  diagram  for  the 
sensitivity  analysis  computer  program  is  shown  in  Appendix  B.  Knowing 
the  critical  DO  in  the  stream,  it  is  then  possible  to  reckon  the  percent 
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deviations  in  the  predicted  DO  from  the  reference  value  when  one  of 
the  reaction  coefficients  takes  different  values  within  its  practical 
range  of  values  while  the  other  reaction  coefficient  assumes  its 
nominal  value  and  all  other  parameters  are  held  at  their  respective 
constant  values* 

Tables  2  to  4  show  the  sensitivity  of  the  hypothetical  river 
system  for  the  variations  in  K2  at  10°C,  20°C,  and  30°C  respectively 
for  the  assumed  initial  conditions  of  BOD  (20  mg/1)  and  DO  deficit 
(1*5  mg/1).  These  are  shown  plotted  in  Figure  2.  When  K2  takes  the 
nominal  value  at  a  particular  temperature,  there  will  be  no  error  in 
predicting  the  response  of  the  system*  However,  when  the  actual  value 
of  K2  is  lesser  or  greater  than  the  nominal  value,  the  actual  value  for 
critical  DO  will  be  different  from  its  nominal  (reference)  value*  The 
figure  clearly  shows  the  temperature  dependence  of  the  sensitivity  of 
the  system  for  the  variations  in  K2*  The  error  in  prediction  increases 
with  increased  deviations  of  the  coefficient  K2  and  this  becomes  in- 
creasingly significant  at  higher  temperatures.  Also  it  should  be  noted 
that  there  are  upper  bounds  for  these  errors  when  the  actual  values  for 
K2  are  lesser  or  greater  than  the  nominal  value*  When  the  value  is  far 
less  than  the  nominal  value*  the  rate  of  deoxygenation  exceeds  the  rate 
of  reaeration  with  the  result  that  the  critical  deficit  could  reach  its 
limiting  value  equal  to  the  DO  saturation  concentration  for  the  given 
temperature*  In  such  a  case*  maximum  deviation  could  only  be  100  per- 
cent, and  this  effect  is  seen  for  the  situation  when  K2  equals  0*2  per 
day  at  30°C.  On  the  other  hand  when  K„  has  a  value  considerably  greater 
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FIGURE  2.  EFFECT  OF  VARIATIONS  IN  K,  VALUES  ON  PERCENT 

ERROR  IN  PREDICTING  THE  CRITICAL  DO  AT  DIFFER- 
ENT TEMPERATURES 
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than  the  assumed  average  value,  the  rate  of  reaeration  will  be  con- 
siderable compared  to  the  deoxygenation  rate  with  the  result  that 
the  critical  deficit  in  the  system  vill  be  the  initial  deficit  it- 
self* This  upper  bound  for  the  percent  error  in  DO  prediction  varies 
with  the  temperature. 

Table  5  shows  the  sensitivity  of  the  system  for  the  vari- 
ations in  Kj  at  10°C,  20°C,  and  30°C.  These  are  plotted  in  Figure  3. 
Figures  U   to  6  show  the  relative  effects  of  the  deviations  in  K^  and 
K£  on  the  sensitivity  of  the  system  at  temperatures  10°C,  20°C,  and 
30°C  respectively.  The  curve  marked  as  "lower  scale"  indicates  the 
sensitivity  of  the  system  when  K2  varies  and  all  other  parameters  of 
the  system  remain  constant.  Likewise  the  curve  marked  as   "upper  scale" 
refers  to  the  sensitivity  of  the  system  when  K,  varies  while  other 
parameters  remain  constant.  It  is  seen  from  these  curves  that  except 
at  10°C,  the  variations  in  the  velocity  coefficients  are  equally  sig- 
nificant in  predicting  the  critical  DO.  At  10°C,  the  error  in  prediction 
when  Kj  takes  values  greater  than  its  nominal  value  is  much  more  than 
the  error  due  to  the  variations  in  K2» 

Generally,  it  is  seen  that  the  errors  in  prediction  due  to 
the  variations  in  the  reaction  coefficients  are  considerable,  with  the 
result  that  there  is  no  justification  for  using  average  values  for 
these  parameters  in  predicting  dissolved  oxygen  concentrations.  Hence 
it  is  postulated  that  by  treating  the  velocity  coefficients  as  variable 
coefficients  Instead  of  as  constant  coefficients,  the  values  predicted 
by  mathematical  modeling  can  be  made  to  approach  the  true  value  for  DO 
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FIGURE  3,  EFFECT  OF  VARIATIONS  IN  K&  VALUES  ON  PERCENT  ERROR 
IN  PREDICTING  CRITICAL  DO  AT  DIFFERENT  TEMPERATURES 
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FIGURE  4.  COMPARISON  OF  THE  EFFECTS  OF  VARIATIONS  IN  K1  AND 
K2  ON  PERCENT  ERROR  IN  PREDICTING  THE  CRITICAL  DO 
AT  10°C 
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FIGURE  5.   COMPARISON  OF  THE  EFFECTS  OF  VARIATIONS  IN  Kt  AND 
K„  ON  PERCENT  ERROR  IN  PREDICTING  THE  CRITICAL  DO 
20°C 
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FIGURE  6.  COMPARISON  OF  THE  EFFECTS  OF  VARIATIONS  IN  Kj  AND 
K,  ON  PERCENT  ERROR  IN  PREDICTING  THE  CRITICAL  DO 
AT  30®C 
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better  than  it  is  possible  now.  The  problem  of  evaluating  the 
assimilative  capacities  of  streams  is  approached  by  simulation 
technique  using  the  Monte  Carlo  method  which  makes  possible  the 
treatment  of  the  parameters  Kj  and  K2  as  variable  coefficients. 
The  concepts  and  procedure  involved  in  these  will  be  developed  in 
subsequent  chapters. 
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IV.  PROBABILISTIC  VARIATIONS  IN 
VELOCITY  COEFFICIENTS 


It  has  been  found  that  in  many  of  the  practical  problems 
which  defy  theoretical  solutions,  an  approximate  but  workable  solution 
could  be  found  using  Monte  Carlo  techniques.  In  these  methods,  it  is 
necessary  to  specify  a  priori  the  probability  distributions  for  the 
random  variations  of  separate  events  which  are  then  merged  into  a 
composite  picture.  Few  attempts  have  been  made  to  define  the  chance 
variations,  in  terms  of  probability  measure,  for  the  initial  BOD  and 
DO  at  the  source  of  pollution.  Montgomery  (1964)  assumed  the  daily 
variations  of  waste  loads,  in  terms  of  BOD,  to  be  normally  distributed 
and  the  assumption  was  based  on  considerations  given  to  the  histories 
from  Gary,  Indiana,  Galesburg,  Illinois,  Dallas,  Texas,  and  Racine, 
Wisconsin.  Whereas,  Thayer  (1966)  in  an  attempt  to  formulate  a  stochastic 
model  for  pollution  and  dissolved  oxygen  in  streams  postulated,  without 
verifying  the  assumption,  that  the  initial  BOD  and  DO  in  a  stream  vary 
independently  according  to  the  binomial  distribution  law.  Though  it  is 
conceptually  feasible,  under  certain  simplified  assumptions,  to  consider 
all  the  separate  events  involved  in  the  mathematical  model  for  stream 
assimilative  capacities  as  random  events,  attention  will  primarily  be 
focused  on  the  variations  in  velocity  coefficients  K^  and  I<2  and  the 
proposed  probabilistic  model  will  be  tested  for  known  initial  conditions. 

VARIATIONS  IN  THE  DEOXYGENATION  COEFFICIENT 

As  has  been  mentioned  earlier,  it  is  well  documented  that  the 
deoxygenation  coefficients  of  domestic  and  industrial  wastes  at  any  given 
location  vary  considerably  and  attempts  have  been  made  to  formulate 
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methods  to  circumvent  the  difficulties  and  inaccuracies  arising  out 
of  these  variabilities  (Churchill  and  Buckingham,  1956;  LeBosquet  and 
Tsivoglou,  1950),   In  order  to  place  the  art  of  predicting  stream 
assimilative  capacities  on  a  more  rational  and  time  tested  basis  than 
those  proposed  by  Churchill  and  Buckingham  (1956),  and  LeBosquet  and 
Tsivoglou  (1950),  it  is  necessary  to  extend  the  concepts  of  Streeter 
and  Phelps  (1925)  taking  variabilities  in  reaction  coefficients  into 
account • 

In  sanitary  engineering  practice,  it  has  been  tacitly  assumed 
that  the  characteristics  of  any  waste,  particularly  those  of  domestic 
origin,  vary  from  time  to  time.  In  most  of  the  controlled  studies,  it 
is  customary  to  use  synthetic  sewage  rather  than  the  waste  generated  in 
an  actual  situation.  One  of  the  parameters  which  define  the  character- 
istics of  a  waste  is  the  rate  at  which  it  is  biologically  oxidized. 
Since  the  characteristics  of  wastes  vary  considerably  from  time  to 
time,  the  rate  parameter  for  the  wastes,  treated  or  untreated,  also 
vary  considerably.  The  general  course  of  deoxygenation  and  the  K^  values 
vary  since  biological  processes  deal  with  the  response  of  living  or- 
ganisms to  their  environment  and  do  not  show  the  same  degree  of  uni- 
formity or  consistency  even  under  standard  conditions,  as  purely 
chemical  processes  do.  Even  though,  this  phenomenon  can  be  attributed 
to  the  viability  and  type  distribution  of  microbial  population  present 
and  the  changes  in  characteristics  of  the  substrate  for  the  microbes, 
it  is  not  possible  to  predict  the  inter-relationship  of  changes  in 
waste  characteristics  and  microbial  response  quantitatively  in  a  deter- 
ministic sense.  Due  to  the  inherent  ambiguities  involved  in  the 
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deoxygenation  process,  where  most  of  the  contributing  factors  are 
subject  to  random  fluctuations,  it  is  very  well  suited  for  being 
treated  under  the  theories  of  probability* 

Ki  Values  of  the  Ohio  River  Samples 

Whenever  one  is  inferring  general  laws  on  the  basis  of 
particular  observations  associated  with  them,  the  conclusions  are 
uncertain  inasmuch  as  the  particular  observations  are  only  a  more  or 
less  representative  sample  from  the  totality  of  all  possible  ob- 
servations. One  way  of  reducing  uncertainty  in  a  solution  is  to  collect 
and  base  it  upon  more  observations.  Data  from  a  1957  survey  of  the 
Ohio  River,  collected  and  published  by  the  United  States  Public  Health 
Service  (1960)  in  connection  with  the  evaluation  of  the  effects  of 
navigational  improvements  on  assimilative  capacity  of  the  Cincinnati 
Fool,  is  one  of  the  most  extensive  and  comprehensive  set  of  all  river 
pollution  survey  data  published.  Extensive  information  on  long-term 
(20°C)  BOD  test  results  for  sewage  treatment  plant  effluent  and  river 
samples  downstream  of  effluent  outfalls,  DO  test  results,  photosyntheti c 
oxygen  production  ("Light  and  Dark  Bottle")  test  results,  etc.,  are 
presented.  For  purposes  of  studying  the  probabilistic  variations  of 
the  rate  coefficient  Kj  in  an  actual  situation,  the  long-term  BOD  data 
on.  river  samples  collected  in  the  reach  between  miles  474.6  and  481.45 
is  used  in  this  study.  All  the  pollutional  loads  emanating  from  Cincin- 
nati enter  the  river  above  this  reach.  The  long-term  BOD  data  on  river 
samples  are  abstracted  from  the  published  data  and  presented  in 
Appendix  D. 
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Having  the  BOD  data  available,  the  next  step  was  to  compute 
the  monomolecular  curves  which  best  fit  the  observed  data*  It  was 
noted  earlier  that,  because  of  its  traditional  usage  and  general  ac- 
ceptance, it  will  be  assumed  that  the  first  stage  of  deoxygenation 
proceeds  according  to  a  monomolecular  or  first  order  reaction*  Some 
method  had  to  be  selected  of  the  numerous  ones  reported  in  the  liter- 
ature for  the  evaluation  of  the  doxygenation  coefficient  and  the  corres- 
ponding ultimate  BOD* 

One  of  the  earlier  methods  for  the  calculation  of  the  velocity 
coefficient  K^  was  that  proposed  by  Reed  and  Theriault  (1931)  which 
uses  the  least  squares  principle  for  minimizing  the  residual  errors 
between  the  observed  data  and  the  fitted  monomolecular  curve*  In  the 
application  of  the  procedure,  a  method  of  successive  approximations 
involving  trial  and  error  is  used  in  which  it  is  necessary  to  assume  a 
trial  value  of  Kj_,  and  this  is  then  used  together  with  the  observed 
data  in  determining  a  correction  factor  to  the  assumed  Kj  value*  The 
operation  is  repeated  with  the  corrected  Kj  value  until  the  correction 
factor  is  small,  or  in  other  words,  until  the  assumed  and  calculated 
values  of  K.  approach  one  another*  Obviously  if  the  initial  estimate 
is  good,  the  number  of  iterations  will  be  minimal,  but  even  if  this 
estimate  is  poor  the  procedure  will  eventually  yield  a  good  estimate 
of  K1#  The  time  necessary  for  hand  computation  is  excessive  and  this 
no  doubt  had  been  a  factor  in  the  development  of  less  rigorous  com- 
putational and  graphical  procedures  (Fair,  1936$  Moore,  Thomas,  and 
Snow,  1950;  Thomas,  1937)*  Since  it  is  generally  considered  (Gannon, 
1963;  Zanoni,  1967)  that  the  Reed -Theriault  method  affords  the  most 
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reliable  estimation  for  K^  and  ultimate  first  stage  BOD,  this  method  is 
used  in  this  study  in  evaluating  these  parameters  for  the  observed  BOD 
data.  The  flow  diagram  for  the  computer  program  is  shown  In  Figure  7, 
and  the  tolerance  limit  between  the  assumed  and  corrected  values  for  K| 
is  set  at  0.005  per  day.  The  terminology  used  in  the  flow  diagram  and 
in  the  computer  program  is  the  same  as  that  adopted  by  Reed  and  Theriault 
(1931). 

In  determining  the  parameters  for  the  first  stage  carbonaceous 
demand,  it  is  necessary  to  make  corrections  for  the  initial  lag  in 
oxygen  uptake  and  oxygen  demand  due  to  nitrification,  if  these  effects 
are  significant  in  the  BOD  progression  data.  Figure  8  shows  the  long- 
term  (20°C)  BOD  data  for  effluent  samples  from  Mill  Creek  and  Little 
Miami  sewage  treatment  plant.  It  is  readily  seen  that  there  is  no  lag 
in  initial  oxygen  uptake  in  either  of  these  samples.  Also,  the  trend 
in  BOD  progression  deviates  on  or  around  the  tenth  day,  heralding  the 
onset  of  nitrification  in  the  bottle  tests.  Figure  9  shows  the  long- 
term  (20°C)  BOD  progression  results  for  the  river  samples  obtained  in 
the  Ohio  River-Cincinnati  Pool  area.  Here  again,  no  noticeable  lag  is 
discerned  and  the  nitrification  process  is  found  to  set  in  on  or  around 
the  tenth  day.  The  BOD  progression  curves  presented  here  are  typical 
of  the  results  obtained  in  the  Ohio  River  survey.  Since  the  flow-through 
time  involved  in  the  survey  is  never  more  than  four  days,  oxygen  demand 
due  to  nitrification  is  not  considered  in  this  study.  Also,  the  BOD 
progression  data  observed  up  to  the  tenth  day  were  used  in  the  Reed- 
Theriault  method  for  computing  Kj  and  ultimate  first  stage  BOD  values. 
In  most  cases  there  were  six  observations  within  ten  days  in  each  series 
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Read  input  values  for  BOD  (Y),  time  (t) 
days,  estimated  value  for  K^y  and  tolerance 
limit  for  correction  factor. 


Write  output 
values  of  L 
and  K, 


I  -  1 


f^  - 


Compute 
[l  -  exp  (-Kjt)] 
t  exp  (-Kjt) 


Ef2' 


"iT 


EfiLi 


r£l£2. 


zf?v 


If2^   -   (Ef1f2)' 


H 


Correction  factor,   H 
Zfx2   E  f2Y  -    Efjf, 


EfiY. 


Ef. 


EfoZ  -  (E  f 


"7TT-J 


Ki^L 


Write  error  message   « 


Yes 


End 


FIGURE   7.     FLOW  DIAGRAM  FOR  COMPUTING  K,   AND  ULTIMATE  FIRST  STAGE  BOD 

VALUES   USING  REED-THERIAULT  METHOD 
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of  BOD  progression  tests  and  five  observations  in  a  very  few  cases. 
The  values  for  Kj  (20°C)  and  ultimate  first  stage  BOD  computed  for 
each  series  of  observations  on  the  Ohio  River  samples  are  shown  in 
Appendix  D. 

Randomness  Test 

Many  nonparametric  statistical  procedures  are  based  on  the 
assumption  of  a  random  sample  of  univariate  observations.  Also 
randomness  of  a  set  of  observations  is  sometimes  itself  the  property 
investigated.  Although  the  concept  of  independent  observations,  a 
random  sample,  is  well-known,  it  is  worthwhile  to  define  it  here.  A 
sample  of  n  observations,  Xj,  %2» •••»xn   from  a  population  with  dis- 
tribution function  F(x)  is  called  a  random  sample  from  that  population 
if 

P  {x1,x2,...,xn}  -  P  {Xl}  P  {x2}   ...  P  {xj       (13) 

where  P  {•}  is  the  probability  of  observing  a  value  less  than  or  equal 
to  the  stated  value.  When  a  set  of  observations  is  analysed,  the  first 
step  is  to  investigate  whether  the  observations  may  be  regarded  as  a 
random  sample  or  not.  The  next  step  is  to  formulate  a  hypothesis  for 
the  type  of  distribution  of  the  population  from  which  the  sample  is 
supposed  to  be  drawn  randomly  and  to  estimate  values  of  the  parameters 
of  the  assumed  distribution  based  on  the  sample. 

The  most  common  case  attempted  is  that  showing  the  relationship 
of  an  effect  to  many  causes  of  which  a  small  number  of  the  causes  exert 
greater  influence  than  do  all  others.  In  such  a  case,  when  neglected 
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vari Abies,  inherent  errors,  and  nonhomogeneity  of  data  have  relatively 
small  effects,  the  relationship  between  the  remaining  limited  number 
of  variables  would  indicate  a  narrow  spread  around  a  basic  function, 
Stoltenberg  and  Sobel  (1965)  in  their  study  to  determine  the  effects 
of  sampling  location,  organic  loading,  tidal  stage,  and  sample  col- 
lection temperature  on  K^  values  for  various  samples  obtained  in  the 
Delaware  Estuary,  came  to  the  conclusion  based  on  analysis  of  covari- 
ance  that  these  factors  did  not  have  significant  effect  on  K^  values. 
Hence  it  is  hypothesized  that  the  variations  in  K^  values  are  generally 
random,  not  attributable  to  any  definite  cause  with  certainty  and 
this  conjecture  was  subjected  to  statistical  randomness  tests  using 
the  observed  K.  values  for  the  Ohio  River  samples.   If  the  hypothesis 
fails,  then  one  has  to  separate  the  trend  in  the  series  of  observations 
made  and  study  the  random  variations  about  the  observed  trend. 

Determination  of  Sequence  Order.  The  first  step  in  applying 
a  test  of  the  random  sample  hypothesis  consists  of  placing  the  n  ob- 
servations considered  in  an  appropriate  sequence  order  (Walsh,  1962). 
The  method  of  assigning  the  integers  l,...,n  to  the  n  observations 
depends  on  the  experimental  situation  considered.  A  common  method  of 
assignment  is  on  the  basis  of  the  times  at  which  the  observations  are 
produced.  Then  the  first  observation  produced  is  assigned  the  integer 
1,  etc.  The  observations  can  be  ordered  on  the  basis  of  location. 
Also  a  combination  of  location  and  time  can  be  used.  In  the  Ohio 
River  survey,  samples  from  the  river  reach  under  consideration  were 
collected  for  analysis  from  five  locations  starting  from  the  head  end 
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of  the  reach  and  moving  downstream  to  cover  all  the  locations  and  the 
whole  operation  was  repeated  several  times.  Since  the  fate  of  pol- 
lution and  the  condition  of  oxygen  concentration  in  the  river  as  the 
waste  discharge  flows  along  the  river  are  of  primary  importance,  the 
observations  for  K.  values  are  ordered  on  the  basis  of  the  times  at 
which  the  observations  are  produced.  The  sequence  of  observations 
thus  obtained  was  then  subjected  to  randomness  tests  which  are  dis- 
cussed in  the  following  sections. 

Another  approach  would  be  to  group  the  observations  according 
to  the  place  of  observation  so  that  any  decreasing  trend  in  the  Kj 
values  could  be  discerned  and  isolated.  It  is  likely  that  the  Kj  values 
of  river  samples  decrease  in  the  downstream  sections  as  the  more 
readily  degradable  organics  are  oxidized.  This  approach  is  not  con- 
sidered in  this  study  since  the  river  sections  from  which  samples  were 
obtained  in  the  Ohio  liver  (USPHS,  1960)  are  not  separated  far  enough 
to  give  any  meaningful  result. 

Runs  Up  and  Down  Test.  Considering  a  sequence  of  n  different 
observations  x^,  X2»...,xn  and  the  sequence  of  signs  (+  or  -)  of  the 
(n  -  1)  differences  (x...  -  x. ),  a  sequence  of  successive  +  signs  is 
called  a  run  up  and  a  sequence  of  successive  -  signs  a  run  down  (Hald, 
1960;  Walsh,  1962).  The  length  of  a  run  is  given  by  the  number  of  same 
signs  defining  the  run.  The  total  number  of  runs  is  denoted  by  R,  the 
number  of  runs  of  length  i  by  r^  and  the  number  of  runs  of  length  k  or 
more  by  Rk  where 


n  -  1 

Z 
i  -  k 


Rk  -    Z  ri  (14) 
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The  expected  value  for  the  total  number  of  runs  E(R),  and  the  variance 
In  the  number  of  runs  V(R)  are  given  by 

E(R)  -  i  (2n  -  1)  (15) 

and 


V(R)  -  ^  (16n  -  29)  (16) 


Also  for  n  >  20,  R  may  be  regarded  as  normally  distributed 
with  good  approximation*  The  expected  number  of  runs  up  and  down  of 
length  k  or  more  in  random  arrangements  of  n  different  observations  for 
values  of  k  from  1  to  7  are  given  by  Hald  (1960)  in  Table  13.10. 

If  the  hypothesis  alternative  to  the  hypothesis  of  randomness 
is  that  a  gradual  change  in  the  level  of  distribution  has  taken  place 
during  the  drawing  of  the  samples,  such  a  change  will  produce  trends 
or  cycles  in  the  observations,  so  that  one  or  more  long  runs  may  be 
expected  to  occur  and  the  total  number  of  runs  will  be  small.  The 
hypothesis  of  randomness  may  therefore  be  tested  by  means  of  the  total 
number  of  runs,  a  small  number  of  runs  being  significant  and  further 
by  means  of  length  of  runs,  very  long  runs  being  significant. 

The  sequence  of  observations  for  Kj^  values  gives  the  dis- 
tribution of  runs  up  and  down  as  shown  in  Table  6.  For  a  total  number 
of  observations  of  83  for  K,  values,  the  expected  number  of  runs  and 
the  corresponding  variance  are  given  by  equations  15  and  16  as  55  and 
14,4  respectively.  The  difference  between  the  observed  and  expected 
number  of  runs  is  -2.  The  standard  deviation  for  the  total  number  of 
runs  is  3.8,  the  variance  being  14.4.  Hence  the  deviation  of  the 
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TABLE  6 

DISTRIBUTION  OF  RUNS  UP  AND  DOWN  FOR  THE 

SEQUENCE  OF  OBSERVATIONS  FOR  Kx  VALUES 

OF  THE  OHIO  RIVER  SAMPLES 


Length 
of  Run 
i 

Number 
Up 

of  Runs  Observed 

Both 
Down    (x-j)    Ry. 

Expected 
E(r|) 

Number 
E(Rk) 

1 

18 

15 

33 

53 

34.7 

55 

2 

5 

8 

13 

20 

16.3 

20.3 

3 

3 

2 

5 

7 

2.9 

4.0 

4 

0 

2 

2 

2 

0.9 

1.1 

>  5 

0 

0 

0 

0 

0.2 

0.2 

Total 

26 

27 

53 

■BIB 

55.0 

*e*» 

observed  total  number  of  runs  from  the  expected  value  is  not  significant. 
Also  very  long  runs  are  not  present  in  the  sequence  of  observations. 
Hence  the  hypothesis  that  the  observations  for  K1  values  form  a  random 
sample  cannot  be  rejected. 

Turning  Points  Test.  This  test  consists  in  counting  the 
number  of  peaks  and  troughs  in  the  series  (Kendall  and  Stuart,  1966), 
A  "peak"  is  a  value  which  is  greater  than  the  two  neighboring  values. 
If  there  are  two  or  more  equal  values  which  are  greater  than  their 
predecessor  and  successor  (a  rare  event  in  general)  they  will  be  regarded 
as  defining  one  peak.  Likewise  a  "trough*4  is  a  value  which  is  lower 
than  its  neighbors.  Both  peaks  and  troughs  are  treated  as  cases  of 
"turning  points"  and  the  interval  between  two  turning  points  is  called 
a  phase.  The  expected  number  of  turning  points,  P,  variance  inP,  V(j»), 
and  the  expected  number  of  phases  Pj  of  length  i  for  a  sequence  of  n 


>bservations  are  given  respectively  as 
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E(P)  -  2  (n  -  2) 


(17) 


V(P)  -  ^  (16n  -  29) 


(18) 


ind 


E(P  )     -  2(n  -  ^-  2)(^2  +  3j  +  1) 
MPi;  (i  +  3)J 


(19) 


For  the  river  K.  values,  the  actual  number  of  turning  points 
in  the  sequence  of  observations  is  52  as  against  the  expected  value  of 
>4,7.  The  distribution  of  "phases'*  for  the  observed  Kj  values  are  shown 
In  Table  7.  The  observed  values  for  the  total  number  of  turning  points 
ind  for  the  "phases"  of  different  lengths  are  so  close  to  the  theoretical 
ralues  that  the  hypothesis  of  random  variations  in  K.  values  cannot  be 


rejected. 


TABLE  7 

DISTRIBUTION  OF  "PHASES"  FOR  THE  SEQUENCE 
OF  OBSERVATIONS  FOR  Kj  VALUES 
OF  THE  OHIO  RIVER  SAMPLES 


Phase 

Number 

of 

Phases 

Number  of  Phases 

Length 

Observed 

Theoretical 

1 

32 

33.3 

2 

13 

14.5 

3 

5 

4.1 

4 

1 

0.9 

Total 

51 

52.8 
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Hypothesis  and  Hypothesis  Testing 

The  frequency  of  occurrence  of  the  observed  K^   (20°C)  values 
for  the  Ohio  River  samples  in  different  class  intervals  are  shown  in 
Table  8,     The  sample  mean  and  standard  deviation  are  respectively  0,173 
per  day  and  0.066  per  day  with  a  coefficient  of  variation  of  38.0  per- 
cent.    The  histogram  for  the  observed  frequency  distributions  is  shown 
in  Figure  10  along  with  a  theoretical  normal  curve  having  a  mean  and 
standard  deviation  of  0.173  per  day  and  0.066  per  day.     Figure  11  shows 
these  results  plotted  on  a  normal  probability  paper.     It  is  seen  that 
all  the  observed  points  lie  very  close  to  a  straight   line  except  at  the 
tail  ends  of  the  curve  fitting  where  the  deviations  of  the  points  from 
the  straight  line  are  negligibly  small*     Since  the  observed  frequency  for 
K^  values  fits  very  well  with  the  theoretical  normal  distribution,   it  is 
hypothesized  that  K^  values  for  Ohio  River  samples  at  Cincinnati  Pool, 
vary  randomly  according  to  Gaussian  probability  with  a  mean  of  0.173  per 
day  and  a  standard  deviation  of  0.066  per  day.     This  hypothesis   is  sub- 
jected to  further  statistical  tests  as  outlined  in  the  following  sections, 

TABLE  8 

FREQUENCY  DISTRIBUTION  OF  OBSERVED  KL  VALUES 
OF  THE   OHIO  RIVER  SAMPLES 


Class  Interval 

Number 

Percent 

Cumula- 

of Ob- 

Relative 

Relative 

tive 

servations 

Frequency 

Frequency 

Frequency 

0.000  to  0.049  1 

0.050  to  0.099  11 

0.100  to  0.149  20 

0.150  to  0.199  23 

0.200  to  0.249  15 

0.250  to  0.299  11 

0.300  to  0.349  2 


0.012 

1.2 

1.2 

0.133 

13.3 

14.5 

0.241 

24.1 

38.6 

0.277 

27.7 

66.3 

0.181 

18.1 

84.4 

0.133 

13.3 

97.7 

0.023 

2.3 

100.0 
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FIGURE  10.     FREQUENCY  PLOTS  FOR  Kt  VALUES  OF  THE  OHIO  RIVER  SAMPLES 
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FIGURE  11.     CUMULATIVE  PROBABILITY  DISTRIBUTION  OF  K,  VALUES 
FOR  THE  OHIO  RIVER  SAMPLES 
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Chi -Square  Test.  This  is  a  simple  test  generally  used  for 
testing  whether  the  data  constitute  a  sample  from  a  population  with 
probability  density  function  (p.d.f .)  f (x)  at  the  significance  level  a 
(Ostle,  1964).  The  range  of  observed  values  is  divided  into  a  number 
of  categories;  and  the  expected  number  of  samples  ej,  that  will  fall 
in  each  category  i,  as  predicted  by  the  p.d.f.  f(x),  is  calculated. 
The  observed  frequencies  nj  in  each  of  the  categories  are  compared 
with  the  expected  frequencies  by  computing  the  chi -square  statistic 
defined  by 

r   (n,  -  e.)2 

X2  «    E     h (2°) 

i  -  1 

Let  f  be  the  number  of  degrees  of  freedom,  given  by  f  - 
r  -  1  -  g  in  which  r  is  the  total  number  of  categories,  and  g  is  the 
number  of  population  parameters  estimated  from  the  sample  and  used  in 
theoretical  distribution.  If  the  calculated  values  of  X2  exceeds  the 
value  of  x  a  f»  the  hypothesis  that  the  sample  belongs  to  the  distribution 
f(x)  is  rejected  at  significance  level  a. 

The  result  of  chi -square  test  is  shown  in  Table  9.  Merging 
the  first  and  the  last  of  the  categories  indicated  in  the  table  with 
their  respective  neighboring  categories,  the  resulting  number  of  cate- 
gories is  5  and  the  number  of  degrees  of  freedom  is  3.  Since  x2  "  1*878 
being  less  than  the  critical  value  x25  3  -  9,35,  the  hypothesis  cannot 
be  rejected. 


TABLE  9 


CHI -SQUARE  TEST  FOR  THE  GOODNESS  OF  FIT 
CONCERNING  THE  HYPOTHESIS  FOR  OBSERVED  Kl   VALUES 
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Class  Interval   Observed   Expected 


0.000  to  0.049 
0.050  to  0.099 
0.100  to  0.149 
0.150  to  0.199 
0.200  to  0.249 
0.250  to  0.299 
0.300  to  0.349 


Frequency  Frequency  n,  -  e^   (nj  -  e.)  /ej 
(nj)       (ej) 


1 

2.6 

11 

8.5 

0.9 

0.073 

20 

19.0 

1.0 

0.053 

23 

24.6 

1.6 

0.104 

15 

18.3 

3.7 

0.748 

11 

7.8 

3.0 

0.900 

2 

2.2 

2  -  1.878 


Kolmogorov-Smirnov  Test.     An  alternative  to  the  chi -square 
goodness  of  fit  test  is  provided  by  the  Kolmogorov-Smirnov  test   (Massey, 
1951;   Ostle,   1964).     If  a  population  is  thought   to  have  some  specified 
cumulative  distribution  function  (c.d.f .)  F(x),   that   is,   for  any 
specified  value  of  x,   the  value  of  F(x)   is  the  proportion  of  individuals 
in  the  population  having  measurements   less  than  or  equal  to  x.     The 
cumulative  step-function  of  a  random  sample  of  n  observations  is  ex- 
pected to  be  fairly  close  to  this  specified  distribution  function.     If 
it  is  not  close  enough,    it  is  evident  that  the  hypothetical  distribution 
is  not  the  correct  one. 

If  Sn(x)   is  the  observed  cumulative  step-function  of  a  sample, 
I.e.,   sn(x)  -»  k/n,   where  k  is  the  number  of  observations   less  than  or 
equal  to  x,   then  the  sampling  distribution  of  d  being  the  maximum  of 
IF(x)   -  Sn(x)  I  is  taken  as  the  statistic  for  testing  the  goodness  of 
fit.     The  d  statistic  should  be  less  than  the  critical  value  which  de- 
pends on  the  number  of  observations  on  which  the  hypothesis  concerning 
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the  population  distribution  is  based.  Massey  (1951)  has  presented 
proof  for  the  fact  that  the  Kolmogorov-Smirnov  test  is  more  powerful 
than  the  chi -square  test  and  that  this  test  will  detect  smaller  devia- 
tions in  cumulative  distributions  than  the  chi-square  test. 

Results  of  the  d  test  for  the  Kj  values  observed  on  the  Ohio 
River  samples  are  shown  in  Table  10.  From  the  table  it  is  seen  that 
the  maximum  value  for  the  difference  between  the  theoretical  and  observed 
cumulative  distributions  is  0.035.  The  critical  value  for  the  d  sta- 
tistic at  5  percent  significance  level  for  a  total  of  83  observations 
is  0.149.  Since  the  observed  value  for  the  d  statistic  is  less  than 
the  critical  value,  the  hypothesis  cannot  be  rejected.  The  Kolmogorov- 
Smirnov  test,  applied  to  individual  observations  instead  of  the  grouped 
data  as  presented  in  Table  10,  is  also  satisfied. 

TABLE  10 

KOLMOGOROV-SMIRNOV  TEST  FOR  THE  GOODNESS 

OF  FIT  CONCERNING  THE  HYPOTHESIS 

FOR  OBSERVED  Kx  VALUES 


Value  of 
Variable 

the 
(x) 

Sn(x)  - 

k 
n 

Normalized 
Variable 

F(x) 

IF( 

x)    -  Sn(x)l 

0.05 

0.012 

-1.864 

0.031 

0.019 

0.10 

0.145 

-1.106 

0.134 

0.011 

0.15 

0.386 

-0.348 

0.368 

0.023 

0.20 

0.663 

+0.409 

0.659 

0.004 

0.25 

0.844 

+1.167 

0.879 

0.035 

0.30 

0.977 

+1.924 

0.973 

0.004 

0.35 

1.000 

+2.682 

0.996 

0.004 

Since  the  effect  of  temperature  on  the  deoxygenation  coefficient 
K}  has  been  extensively  studied  and  found  to  be  adequately  expressed  by 
the  deterministic  equation,  Eq.  8,  it  is  assumed  that  the  variations  in 
Kj  values  at  any  temperature  other  than  20°C  will  also  be  normally 
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distributed.  Also  no  attempt  was  made  in  this  study  to  verify  the 
goodness  of  fit  with  other  theoretical  distributions  like  lognormal 
distribution,  gamma  distribution,  t  and  F  distributions,  etc.  However, 
the  procedure  involved  in  the  verification  of  the  goodness  of  fit  with 
any  other  theoretical  distribution  will  be  the  same  as  that  for  the 
normal  distribution  discussed  in  this  work. 

VARIATIONS  IN  REAERATION  COEFFICIENT 
Background 

The  rate  of  reaeration,  under  constant  conditions  of  tempera- 
ture and  turbulence,  is  directly  proportional  to  the  oxygen  deficit  in 
the  water.  For  many  years  attempts  have  been  made  to  derive  a  method 
which  would  permit  the  prediction  of  reaeration  coefficient,  relating 
some  of  the  easily  obtainable  hydraulic  characteristics  of  a  stream. 

O'Connor  and  Dobbins  (1956)  presented  a  theoretical  derivation 
of  the  reaeration  coefficient  in  which  they  attempted  to  define  the  rate 
of  reaeration  in  terms  of  the  rate  of  renewal  of  the  surface  film,  and 
assumed  that  the  best  estimate  was  given  by  the  ratio  of  the  vertical 
velocity  fluctuation  and  mixing  length.   In  streams  in  which  there  is 
pronounced  velocity  gradient,  the  following  equation  for  nonisotropic 
turbulence  was  developed: 

480  D„0-5  S0-25 

*2  "  — J5 

and  for  isotropic  turbulence,  In  relatively  deep  channels  where  there 
is  no  pronounced  velocity  gradient: 
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127  D  °'5  S0'25 
K2  "  fl (22) 


where  Dm  is  the  coefficient  of  molecular  diffusion,  S  is  the  channel 
slope  and  H  is  the  mean  depth.  Isotropic  turbulence  was  assumed  ar- 
bitrarily for  Chezy  coefficients  greater  than  17  and  nonisotropic 
turbulence  for  those  less  than  17.  O'Connor  and  Dobbins  (1956)  at- 
tempted to  verify  their  equations  by  comparing  predicted  values  of  the 
reaeration  coefficients  for  the  Ohio  River  using  the  Streeter-P helps 
oxygen  sag  equation  (1925)*  They  obtained  good  agreement  in  spite  of 
the  fact  that  errors  in  the  estimation  of  the  parameters  such  as  the 
time  of  water  travel  through  each  reach,  the  DO  concentration,  ultimate 
BOD,  deoxygenation  coefficient,  etc.,  for  use  in  the  Streeter-Phelps 
equation,  are  reflected  in  the  computed  value  of  K£. 

Churchill  et  al.  (1962),  using  multiple  regression  techniques, 
determined  the  coefficients  of  different  equations,  formulated  by  dimen- 
sional analysis,  relating  observed  reaeration  rates  and  various  hydraulic 
parameters  of  river  reaches.  Corrections  for  the  effects  of  photosyn- 
thesis wen  made  in  computing  reaearation  coefficients  from  observed 
measurements  of  changes  in  oxygen  concentration  in  selected  stream 
reaches.  They  concluded  that  it  is  possible  to  predict  reaeration  rate 
for  a  river  reach  from  the  average  velocity,  depth,  and  temperature  of 
the  water  using  the  equation: 


K2(20)  "  5'026  V0'969  H"1'673  (23) 
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where  V  is  the  mean  velocity  of  flow  and  H  is  the  mean  depth.  They 
also  suggest  that  suitable  correction  should  be  made  for  the  effects 
of  pollution  in  predicting  stream  K2  values.  Inclusion  of  other  stream 
characteristics  such  as  slope,  friction  coefficient,  and  Reynolds 
number  did  not  markedly  increase  the  accuracy  of  the  predicted  value. 
A  direct  correlation  between  experimentally  determined  re- 
aeration  coefficients  in  an  artificial  channel  and  the  corresponding 
values  for  longitudinal  mixing  coefficient  was  demonstrated  by  Krenkel 
and  Orlob  (1963).  They  proposed  the  empirically  derived  equation: 

K2(20)  "  1«138  x  10~5  Dl1,321  H"2'32  (24) 

where  EL  is  the  longitudinal  mixing  coefficient  and  H  is  the  depth  of 
flow. 

Dobbins  (1964)  proposed  a  mathematical  model  for  the  re- 
aeration  process  based  on  the  concept  of  "an  interfacial  liquid  film" 
which  maintains  its  existence  in  the  statistical  sense,  that  is  the 
film  is  always  present  but  the  liquid  content  of  the  film  is  being 
continuously  replaced  in  a  random  manner  by  the  liquid  from  the  main 
body.  This  concept  led  to  the  formulation  of  the  equation 


K 


i       rLt 
-  (D^r)*  coth  qgjp  (25) 


in  which  KL  is  the  liquid  film  coefficient,  r  is  the  average  frequency 
with  which  the  stagnant  film  is  replaced  and  Lt  is  the  thickness  of 
film.  He  further  developed  several  equations  using  postulated  constants 
to  arrive  at  the  relationship  between  K^,  r,  L,  and  K2  applicable  to 
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natural  streams.  His  attempts  to  establish  numerical  values  for  the 
hypothesized  constants  and  proportionality  factors  were  somewhat  less 
than  successful.  Thackston  and  Krenkel  (1965)  have  raised  objection 
to  the  correctness  of  the  theory  proposed  by  Dobbins  (1964),  questioning 
the  validity  of  the  propounded  relationships  between  the  constants  used 
in  the  theory,  some  of  which  are  based  solely  on  judgment. 

Several  other  attempts  to  formulate  empirical  relationships 
for  the  prediction  of  K2  values  in  rivers  have  been  reported  in  the 
literature.  Owens  et  al.  (1964)  proposed  an  equation  based  on  the  data 
efatatsed  in  several  reaches  of  lowland  and  Lake  District  streams  in 
England,  using  multiple  regression  analysis  as: 

K2(20)  "  10«09  v°*73  H"1'75  (26) 

Also  Langbein  and  Durum  (1967)  gave  the  following  relationship: 

K2  -  3.3  V  H"1*33  (27) 

In  all  these  prediction  equations  developed  either  on  theoret- 
ical basis  or  using  regression  analysis,  considerable  scattering  of  the 
observed  values  for  K2  are  noticed  about  the  prediction  line.  As  has 
been  mentioned  earlier,  in  one  of  the  reaches  of  the  Holston  River  in 
the  Tennessee  Valley,  the  observed  K2  value  ranged  from  0.10  to  1.18  per 
day  at  20°C  for  constant  discharge  in  the  river  whereas  the  predicted 
value  for  the  reach  was  0.59  per  day*  Also  in  another  reach  of  Watauga 
River  which  is  a  shallow  and  rapid  stream  compared  to  Holston  River, 
observed  values  for  K2  varied  from  2.26  to  8.86  per  day  whereas  the  pre- 
dicted value  using  regression  Eq.  23  developed  for  these  rivers  was  3.22 
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per  day.  Most  of  the  mathemat leal  models  developed  so  far  for  predicting 
K2  values  in  rivers  use  at  best  the  average  values  for  river  velocity 
and  depth.  Since  these  change  from  section  to  section  and  also  within 
the  section  itself  in  natural  streams,  the  estimated  average  values  do 
not  yield  satisfactory  results.  Also  errors  in  sampling  and  analyses, 
and  in  the  estimation  of  mean  travel  time  combined  with  the  uncertainties 
of  wind  effects,  all  complicate  the  prospects  of  predicting  K2  reasonably 
accurately. 

Eckenf elder  and  O'Connor  (1961)  indicated  that  the  variations 
in  K2  depend  upon  the  hydraulic  properties  of  the  particular  channel: 
roughness,  width  and  curvature,  and  that  the  variations  are  usually 
defined  by  a  normal  distribution.  They  further  suggested  that  any  value 
of  the  reaeration  coefficient  within  limited  statistical  ranges  could  be 
used  in  oxygen  balance  calculations.  O'Connor  (1958)  presented  data, 
taken  from  actual  stream  surveys,  showing  variations  of  depths  from 
station  to  station  in  the  Wabash,  Clarion  and  Codorus  Rivers.  These 
were  found  to  be  normally  distributed  and  he  discussed  the  importance 
of  considering  these  variations  in  determining  the  required  number  of 
cross  sectional  areas  to  insure  for  a  given  probability  that  the  measured 
mean  depth  for  use  in  prediction  equations  for  K2  will  fall  within  a 
given  percent  of  the  true  mean. 

Not  only  the  variations  in  depth  affect  the  evaluation  of  K2 
values,  but  also  several  other  factors  discussed  earlier  come  into  play. 
It  becomes  necessary  to  ascertain  the  combined  effects  of  all  these 
factors,  since  it  is  not  possible  to  assess  the  effects  individually  in 
most  cases.  The  data  collected  by  the  Tennessee  Valley  Authority  (1962) 
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on  natural  streams  with  constant  discharges  for  each  set  of  observations, 
being  the  most  extensive,  will  be  used  in  this  study  to  determine  the 
extent  of  variations  and  the  parameters  characterizing  these  variations. 

Modified  Regression  Equation 

The  data,  obtained  in  Tennessee  Valley  rivers  (1962),  com- 
prised of  509  Individual  observations  under  30  different  sets  of 
experimental  conditions.   In  deriving  the  regression  equations, 
Churchill  et  al.  (1962)  used  geometric  means  of  observations  for 
K2  in  each  set  of  experiments.  In  the  knowledge  of  this  author,  most 
of  the  statistical  theories  concern  themselves  with  the  arithmetic 
means  and  not  with  geometric  means.  In  the  application  of  statistical 
and  probability  theories  to  real  phenomena  two  results  play  a  con- 
spicuous role  (Ostle,  1964;  Parzen,  1966).  These  results  are  known 
as  the  Law  of  Large  Numbers  and  Central  Limit  Theorem,  and  deal  with 
the  arithmetic  means  of  observations.  Hence  a  regression  analysis  of 
the  data  obtained  in  Tennessee  Valley  rivers  using  arithmetic  means 
of  group  observations  was  carried  out  in  this  study.  An  abstract  of 
published  data,  arithmetic  means  of  observations,  and  the  residual 
errors  in  the  regression  analysis  are  shown  in  Appendix  E. 

A  general  linear  regression  equation  may  be  represented  by 

Zq     -  AQ+  klZl   +  AjZj*  ...  ApZp  (28) 

where  Zq   is  the  dependent  variable  estimated  from  the  independent 
variables  Z^,   Z2,...,Z  .  From  the  historical  data,  Z4,  with  i  -  l,...,p 
are  known  and  ZQ*f   the  actual  value  of  dependent  variable  Z0  is  also 
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known.  For  the  least  squares  linear  regression,  the  coefficients  Aq, 
Ap...,A  are  so  evaluated  that  the  sum  of  the  squared  errors 


SE2  -  £  (Z0«  -  Z)2  (29) 


is  minimized.  Details  of  solution  are  available  from  standard  text 
books  on  statistics  (Ezekiel  and  Fox,  1959;  Ostle,  1964).  Z0'  -  Zq 
is  referred  to  as  the  residual  error,  or  as  the  random  component, 
and  so  the  least  square  regression  minimises  the  variance  of  the 
random  component. 

The  coefficients  A^,  with  i  -  l,...,p  are  known  as  the  net 
multiple  regression  coefficients  of  Zq  on  Zj,  The  coefficient  of 
multiple  correlation  Rq.1  2       is  a  measure  of  tne  proportion  of 
the  variance  of  Zq  that  is  explained  by  the  multiple  linear  regression 
equation.  It  is  given  by 

P  n 

o  i   *,   Aif,  2  1AZ0JAZ1J] 

R2  -  i * 3 * (30) 

u.i,i,...,P  n    o 

Z      LZ   0j 
J  "  1 

In  the  IBM  7094  system  available  at  the  University  of  Illinois,  there 
is  a  standard  program  of  the  Statistical  Services  Research  Unit  for 
calculating  various  statistical  estimates  of  multiple  linear  regression 
including  the  regression  coefficients,  the  coefficient  of  multiple 
correlation,  etc.  It  was  used  in  this  study  for  multiple  regression 
analysis  of  the  data  for  river  reaeration  coefficients. 
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The  regression  analysis  yielded  an  equation 

K2(20)  -  5.827  V0'924  H-1'705  (31) 

with  a  correlation  coefficient  of  0.917  which  is  better  than  the 
reported  correlation  coefficient  of  0.822  (Churchill  et  al.,  1962) 
using  geometric  means  for  the  group  data.  The  form  of  the  regression 
equation  was  linearized  by  taking  logarithm  of  both  the  sides  of  the 
equation  before  applying  the  general  linear  equation,  Eq.  28.  Results 
of  the  regression  analysis  are  shown  in  Figure  12.  There  is  consider- 
able deviation  of  the  observed  values  for  K2  from  the  predicted  values, 

An  analysis  of  the  distributions  of  percent  error  indicated 
that  they  are  normally  distributed  with  mean  zero  and  a  standard 
deviation  of  36.8  percent  as  verified  by  Kolmogorov-Smirnov  test  for 
goodness  of  fit.  Percent  error  was  obtained  by  dividing  the  residual 
error  by  the  predicted  value  and  multiplying  the  resulting  fraction  by 
100.  Figure  13  shows  the  cumulative  probability  distributions  of  per- 
cent residual  errors  in  predicting  K2  values,  plotted  on  a  normal 
probability  paper.  The  nature  of  the  distribution  of  errors  in  the 
prediction  of  K2  values  is  not  affected  by  temperature  changes,  since 
the  basic  assumption  in  any  regression  analysis  is  that  the  errors  in 
prediction  are  normally  distributed. 

Since  the  error  distributions  are  considered  random,  it  is 
hypothesized  that  K2  values  for  Tennessee  Valley  rivers  could  be  pre- 
dicted by  superimposing  a  random  component  governed  by  a  normal  proba- 
bility distribution  with  a  mean  zero  and  a  standard  deviation  of  36.8 
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percent,  onto  a  trend  component  predicted  by  the  regression  equation, 
Eq.  31 •  Since  the  standard  deviation  of  the  percent  error  involved 
In  the  predict!  <ns  of  K2  values  is  considerable,  there  seems  to  be 
little  justification  in  ignoring  this  random  error  component  altogether 
in  evaluating  stream  reaeration  coefficients. 
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V,  MONTE  CARLO  METHOD  AND  SIMULATION  MODEL  FOR 
STREAM  ASSIMILATIVE  CAPACITIES 


MONTE  CARLO  METHOD 
General 

Monte  Carlo  methods  comprise  that  branch  of  experimental 
mathematics  which  is  concerned  with  experiments  on  random  numbers 
(Hammersley  and  Hands comb,  1964;  Hammersley  and  Morton,  1954;  Kahn, 
1956;  Meyer,  1954).  Problems  handled  by  Monte  Carlo  methods  are 
generally  classified  into  two  types;  probabilistic  or  deterministic 
according  to  whether  or  not  they  are  directly  concerned  with  the 
behavior  and  outcome  of  random  processes.  In  the  case  of  probabilistic 
problems  the  simplest  Monte  Carlo  approach  is  to  observe  random  numbers, 
chosen  in  such  a  way  that  they  simulate  physical  random  processes  of 
the  original  problem,  and  to  infer  the  desired  solution  from  the  be- 
havior of  these  random  numbers.  Monte  Carlo  methods  have  found  ex- 
tensive use  in  the  fields  of  nuclear  physics  and  they  have  been  employed 
in  other  fields  of  science  like  chemistry,  engineering,  biology,  medi- 
cine, etc.  The  significance  and  the  concepts  involved  in  Monte  Carlo 
methods  are  brought  out  by  considering  two  examples. 

The  central  problem  in  designing  a  nuclear  reactor  is  to 
determine  the  critical  size  of  the  core.  Except  when  the  geometry  is 
simple  or  the  neutrons  have  small  energy  ranges,  analytical  techniques 
are  difficult  and  recourse  is  often  made  to  Monte  Carlo  techniques.  A 
free  neutron,  put  into  the  core,  performs  a  random  walk  there,  colliding 
with  fissile  nuclei  until  it  either  escapes  from  the  core  or  is  absorbed 
by  a  nucleus.  If  the  former  occurs,  the  neutron  continues  its  random 
walk  in  the  reflector;  if  the  latter  occurs,  a  fission  may  result 
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ejecting  several  free  neutrons,  called  descendants,  each  of  which 
begins  its  own  random  walk.  The  size  of  the  core  is  critical  when 
the  average  number  of  neutrons  in  it  stays  constant*  Three  proba- 
bility distributions  determine  such  a  random  walk  (Hammersley  and 
Morton,  1954): 

1.  The  distance  a  neutron  travels  between  successive 
collisions  is  distributed  exponentially  with  a  mean  called  the  total 
mean  free  path  (varying  from  medium  to  medium)*  At  a  boundary  a 
neutron  continues,  i  n  the  same  direction  with  the  same  energy,  a 
distance  with  a  mean  equal  to  the  total  mean  free  path  in  the  new 
medium* 

2*  There  are  assigned  probabilities  of  the  possible  types 
of  collision,  namely  capture,  scattering,  or  fission  and  their  results 
including  energy  changes* 

3.  The  direction  a  neutron  takes  after  a  collision  is 
often  isotropically  distributed,  but  for  collisions  with  light  nuclei 
is  peaked  in  the  direction  of  the  incident  neutron's  path. 
These  three  distributions  depend  on  the  energy  of  the  neutron  concerned. 
Knowing  them,  one  can  conduct  a  neutron  from  collision  to  collision  by 
sampling,  recording  the  position,  energy,  and  direction  of  its  motion 
and  its  descendants  just  before  each  collision*  A  tally  for  the  number 
of  neutrons  in  the  core  is  kept  due  to  several  parent  neutrons  in  the 
core  as  a  function  of  a  so-called  "census  parameter,"  namely  a  quantity 
measuring  the  current  duration  of  the  process  (e.g.,  the  number  of 
collisions  to  date,  or  the  total  distance  traveled,  etc*)*  This  having 
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been  done  for  various  core  sizes,  one  can  estimate  the  critical  size 
by  interpolating  for  a  constant  tally  between  the  increasing  and  de- 
creasing tallies.  The  performance  of  a  nuclear  reactor  is  simulated 
here  by  choosing  random  numbers  which  represent  the  random  motions  of 
the  neutron  in  it.  In  this  way,  one  can  experiment  with  a  reactor 
without  incurring  the  cost,  in  money,  time,  and  safety,  of  its  actual 
physical  construction. 

As  an  example  for  the  application  of  the  Monte  Carlo  tech- 
niques in  the  field  of  biology,  if  one  wishes  to  study  the  growth  of 
an  insect  population  on  the  basis  of  certain  assumed  vital  statistics 
of  survival  and  reproduction,  a  model  could  be  set  up  with  paper  entries 
for  the  life  histories  of  individual  insects.  To  each  such  individual, 
random  numbers  are  allotted  for  its  age  at  the  births  of  its  offspring 
and  its  death;  and  then  treat  these  and  succeeding  offsprings  likewise. 
Handling  the  random  numbers  to  match  the  vital  statistics,  one  gets 
what  amounts  to  a  random  sample  from  the  population  which  can  be  anal- 
yzed just  as  real  experimental  data  collected  in  the  field. 

Thus  it  is  seen  that  the  Monte  Carlo  method  involves  essen- 
tially manipulations  of  random  numbers  in  the  simulated  model  for 
physical  processes  and  is  based  on  the  assumption  that  in  a  game  of 
chance,  the  expected  outcome  can  be  estimated  in  principle  by  averaging 
the  results  of  a  large  number  of  plays  of  the  game.  Also  methods  em- 
ploying Monte  Carlo  techniques  have  afforded  workable  solutions  in  a 
widely  divergent  area  of  endeavors. 

There  are  essentially  three  steps  involved  in  solving  a 
problem  using  Monte  Carlo  techniques.  They  are: 

1.  Choosing  or  analogizing  the  probability  process. 
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2.  Generating  sample  values  of  the  random  variables* 

3.  Use  of  variance  reducing  techniques  in  the  Monte  Carlo 
predictions. 

The  first  of  these  will  be  considered  in  a  subsequent  section  of  this 
chapter  in  order  to  maintain  the  continuity  of  the  discussion  on  the 
techniques  of  Monte  Carlo  methods.  Attention  will  be  focused  now  on 
the  methods  generally  used  in  developing  sample  values  of  the  random 
variables. 

Random  Sampling  from  Specified  Distributions 

Sources  of  Random  Numbers.  Mechanical  sampling  devices  like 
disks  and  roulette  wheels  have  been  used  for  generating  random  numbers. 
Since  the  wheel  spinning  methods  are  tedious  and  cumbersome,  they  have 
given  way  to  tables  of  random  numbers  and  use  of  computing  machines  for 
developing  a  series  of  random  numbers. 

There  are  several  tables  of  random  numbers  (The  Rand  Corpor- 
ation, 1955;  Tippet,  1927;  Wold,  1954)  which  have  been  obtained  by 
performing  permutations  on  the  high  order  digits  of  mathematical 
tables,  by  compiling  from  lottery  drawings,  etc.,  or  from  some  chance 
devices.  They  have  been  subjected  to  several  statistical  tests  and 
are  found  acceptable  for  general  sampling  use.  Usually  these  tables 
have  uniformly  distributed  decimal  digits,  though  there  are  some  that 
have  Gaussian  or  normal  deviates.  In  such  tables,  the  first  number 
to  be  used  is  generally  arbitrarily  selected,  perhaps  by  opening  a 
page  at  random,  and  with  eyes  closed,  touching  a  number  in  the  table 
with  a  finger.  The  rest  of  the  numbers  are  usually  chosen  in  sequence. 
Random  number  tables  are  very  convenient  for  generating  a  fairly  small 
sequence. 
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Random  Number  Generators.  Several  mechanical  or  electronic 
systems  have  been  developed  for  generating  a  sequence  of  random  numbers 
or  variables  (Tocher,  1963).  The  electronic  machines  are  based  on  the 
principle  of  converting  a  source  of  random  noise  into  a  train  of  pulses 
which  are  used  to  drive  a  cyclic  counter.  The  drive  is  interrupted 
at  fixed  intervals  of  time,  and  the  successive  status  of  the  counter 
gives  the  successive  random  numbers.  Machines  that  use  analog  systems 
in  which  the  random  input  is  self  generated  are  known  as  analog  ran- 
domisers.  A  detailed  discussion  of  the  random  number  generators  is 
given  by  Tocher  (1963). 

Development  of  Pseudo-Random  Numbers.  Though  random  number 
tables  are  available  as  punched  cards,  with  increasing  use  of  digital 
computers,  mathematical  methods  for  generating  *  'pseudo-random"  numbers 
within  the  computing  machines  have  been  developed  in  order  to  eliminate 
the  need  for  extensive  input  of  random  numbers  (Kahn,  1956;  Meyer,  1954; 
Tocher,  1954),  Also  when  it  is  desired  to  repeat  the  calculations  for 
checking  purposes,  the  random  numbers  used  in  the  initial  calculations 
must  be  available  for  the  later  calculations,  and  hence  pseudo-random 
numbers  are  useful  for  generating  a  reproducible  random  number  sequence. 
Though  the  reproducibility  implies  the  possibility  of  prediction  of  the 
sequence,  and  hence  its  nonrandomness,  the  advantages  of  reproducibility 
make  it  desirable  to  investigate  the  possibility  of  generating  approxi- 
mately random  sequences  by  deterministic  processes.  Pseudo-random 
numbers  may  be  defined  as  reproducible  sequence  of  numbers  developed 
by  using  a  deterministic  process  that  behaves  like  a  random  sequence 
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when  subjected  to  certain  standard  statistical  tests.  Methods  for 
generating  pseudo-random  numbers  include  the  mid-square  method,  the 
mid-product  method,  the  congruence  methods  and  other  miscellaneous 
methods*  Only  the  multiplicative  congruence  method  which  is  used  in 
this  study  is  discussed  below  briefly. 

Multiplicative  Congruence  Method, --One  of  the  characteristics 
of  a  pseudo-random  sequence  is  that  the  sequence  is  cyclic  (Meyer, 
1954;  Tocher,  1954).  Thus  attempts  have  been  made  to  generate  cycles 
of  maximum  length  so  that  a  long  nonrepetitive  sequence  is  produced. 
Lehmer  suggested  the  use  of  the  theory  of  numbers  for  developing  such 
long  cyclic  sequences.  The  theory  of  congruences  deals  with  such 
sequences.   In  this  theory,  a  number  Xj,  is  said  to  be  congruent  with 
X2  modulo  M  if  x2  -  Xj  is  divisible  by  M. 

Consider  a  first  order  recursion  equation  of  the  form 

xn  +  1     -  k^  (mod  M)  (32) 

where  k  and  M  are  initially  arbitrarily  chosen.  A  model  given  by 
Eq,  32  is  known  as  a  congruential  multiplicative  model.  Then  for 
this  model, 


x„  »  knx  (mod  M)  (33) 

n       o 


Since,   at  the  end  of  a  cycle  of  length  n, 


x       -     x       -»     kV    (mod  M)  (34) 

n  o  o 


the  value  of  kn  should  satisfy  the  equation 
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kn  S     l  (mod  M)  (35) 

For  use  in  digital  computers  in  which  numbers  are  represented  in  the 
binary  system,  it  is  preferable  to  assume  M  -  2P  where  p  is  an  Integer, 
to  facilitate  division.   It  can  be  shown  (Tocher,  1954)  that  for  M  -  2P, 
the  maximum  possible  cycle  is  of  length  2P  ,  and  the  only  value  of  k 
giving  this  cycle  must  satisfy  the  condition 

k  -  +3  (mod  8)  (36) 


and  x0 


is  to  be  odd.  A  suitable  choice  of  k  is  any  odd  power  of  5  since 
52<1  +  1  .   -3  (mod  8);  with  q  »  1,  2,...        (37) 


Also  p  is  to  be  chosen  as  the  largest  number  of  binary  digits  shown  in 
the  computer.  After  multiplication  by  k,  the  least  significant  p 
digits  of  the  double  length  product  are  used  as  ^  +  j. 

A  subroutine  available  at  the  University  of  Illinois  for 
generating  random  numbers  uses  k  ■  5",  p  -  36,  and  the  initial  value 
of  x0  is  set  to  235  -  1.   It  is  also  possible  to  preset  this  initial 
value  to  any  other  chosen  value  with  the  restriction  that  the  chosen 
number  should  be  odd.  The  random  sequence  when  read  as  a  binary 
fraction  gives  the  uniformly  distributed  random  number  sequence.  The 
subroutine  has  been  tested  for  the  randomness  of  the  sequence. 

Transformation  of  Random  Numbers 

In  many  problems  a  random  sample  from  a  specified  or  known 
distribution,  which  is  generally  different  from  the  uniform  distribution, 
is  required.  A  random  sample  with  any  given  distribution  can  be  obtained 
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(Kahn,  1956;  Tocher,  1954)  from  a  given  sequence  of  uniformly  distributed 
random  numbers  by  using  the  inverse  probability  integral  transformation. 
This  is  based  on  the  basic  fact  that  for  any  distribution,  the  cumulative 
probability  distribution  function  has  a  rectangular  distribution,  uni- 
form over  (0,1).  Let  x  be  the  random  variable  with  a  cumulative  dis- 
tribution function,  denoted  by  c.d.f.,  and  designated  by  F(x).  The 
random  variable  Y  •  F(x)  is  given  by 


x 

Y  -  F(x)  -  J  f(x)dx  (38) 


where  f(x)  is  the  probability  density  function  denoted  by  p.d.f. 

Consider  the  probability  distribution  of  Y.  Let  the  p.d.f. 
and  c.d.f.  of  Y  be  g(y)  and  G(y)  respectively.  Then, 

G(y)  -  Pr  {Y  <  y}  -  Pr  [X  <  x}  -  F(x)   -  Y 

f or  0  <  Y  <  1  (39) 

and 

g(y)  -  d|£^  -  1  for  0  <  Y  <  1  (40) 

dy  —   — 

This  is  the  probability  density  function  for  the  uniform  (0,1)  distribu- 
tion. By  taking  a  uniform  (0,1)  sample  of  Y,  and  taking  the  Inverse 
x  -  F-1(y)  of  the  probability  integral  transformation  Y  -  F(x),  a  sample 
of  X  is  developed.  By  repeating  this  for  as  many  samples  as  required,  a 
sequence  of  random  numbers  having  any  given  distribution  can  be  developed 
from  a  sequence  of  uniform  (0,1)  random  numbers. 
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The  process  of  Inverse  trans format  ion  can  be  done  by  several 
methods  (Tocher,  1954).  A  relatively  simple  method  for  the  inverse 
probability  integral  transformation  is  the  use  of  interpolation  from  a 
table  having  corresponding  values  of  x  and  y,  and  the  use  of  a  suitable 
degree  of  polynomial  for  interpolation.  This  method,  though  unsophis- 
ticated, is  fairly  simple  and  was  used  in  this  study  to  transform  the 
c.d.f.  to  the  random  number. 

Variance  Reducing  Techniques 

The  methods  which  are  most  often  used  to  reduce  variance  in 
Monte  Carlo  problems  are  straightforward  sampling,  systematic  sampling, 
stratified  sampling,  use  of  expected  values,  correlation,  and  Russian 
roulette  and  splitting  (Kahn,  1956;  Meyer,  1954).  The  methods  which 
can  be  used  to  reduce  variance  are  often  dependent  upon  the  probability 
model  and  in  some  cases  on  the  techniques  used  to  generate  values  of 
the  random  variables.  Kahn  (1956)  has  given  an  excellent  discussion 
of  the  general  nature  of  the  techniques  with  their  applications  to 
evaluation  of  integrals.  The  first  of  these  which  has  found  extensive 
use  in  engineering  problems,  and  which  is  used  in  this  study  is  dis- 
cussed below. 

Straightforward  Sampling.  This  technique  is  based  on  the 
premise  that  one  way  or  reducing  uncertainty  in  an  answer  is  to  base 
it  upon  more  observations  (Hammersley  and  Hands comb,  1964;  Kahn,  1956; 
Tocher,  1954).  Broadly  speaking  there  is  a  square  law  relationship 
between  the  error  in  an  answer  and  the  requisite  number  of  samples. 
If  a  sequence  of  n  random  variables  x^,  X2,«..,xn  are  picked  from  the 
P*d.f.  f (x)  and  If  a  random  variable  sn  defined  by  the  equation 
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n 

*n  -  J   2   ■<*!>  <41> 

i  «  1 


and  the  integral 


J  z(x)f(x)dx  (42) 


exist  in  the  ordinary  sense,  zn  will  almost  always  approach  z  as  a  limit 
as  n  approaches  ». 

The  integral  Eq.  42  is  called  the  expected  value  of  the 
function  z(x),  and  z  is  called  an  estimate  of  z.  If  z  ,  the  expected 
value  of  z2(x)  also  exists,  an  estimate  can  be  made  about  the  amount 
that  zn  deviates  from  z  for  large  n.  Denoting  the  variance  of  z(x)  by 
<r,  it  is  given  by  the  equation 


a2    m     (z  -  »)2  (43) 

-  J  (SB  -  Z)2f (X)dx 

-  Jz2f(x)dx  -  2z  j  zf(x)dx 

+  z2  J  f (x)dx 

m      Z2  -  2Z2  +  Z2 

-  z7  -  z2  (43a) 

Applying  the  Central  Limit  Theorem  (Parzen,  1966),  for  large  n 
the  probability  that  the  event  z  -6  <  zn  <  z  +  6  occurs  is  asymptotically 
independent  of  the  exact  nature  of  z(x)  or  f(x)  but  depends  only  on  n 
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and  cr»  It  can  be  deduced  from  the  Central  Limit  Theorem  that 


,  *&  2 

pr  { S„  <  z   +  8}  -  -±j    ?   exp  (-  5_)dx         (64) 


/*  1    "^-t 


and  that  o^s  »  the  variance  of  the  estimator  z_  of  z  is  given  by  the 
n 


equation 


A     -  £  <45> 


zn 


or 


n     /h 


Thus,  it  is  seen  that  there  is  a  square  law  relationship  between  the 
error  in  an  answer  and  the  requisite  number  of  observations.  To 
reduce  the  error  tenfold  calls  for  a  hundred  fold  increase  in  the 
observations,  and  so  on. 

Estimation  of  Sample  Size.   In  any  simulation,  it  is  necessary 
to  know  approximately  the  sample  size  required  before  sampling  is 
started  in  order  that  the  size  is  neither  too  large  to  make  it  very 
expensive  compared  to  the  information  gained,  nor  too  small  to  be 
reliable.  When  the  required  precision  and  confidence  levels  for  the 
reliability  of  prediction  are  given,  the  sample  size  can  be  determined 
by  the  following  approach  (Flagle,  et  al.,  1960;  Massey,  1951). 


85 

Let  3  and  a  percent  be  the  error  and  significance  levels  for 
Che  reliability  of  predicted  values  in  a  simulation  model.  Within  the 
same  confidence  interval,  the  precision  could  be  increased  only  by  in- 
creasing the  sample  size.  Since  the  standard  deviation  and  consequently 
the  tolerance  diminishes  as  the  square  root  of  the  sample  size,  to  in- 
crease the  precision  by  a  factor  of  two,  the  sample  sise  would  have  to 
be  quadrupled.  If  ^fn)   is  the  critical  value  of  the  *d'  statistic  in 
the  Kolmogorov-Smirnov  test  for  goodness  of  fit  at  a  signficiance 
level  of  a,  which  is  defined  as  the  maximum  absolute  difference  between 
sample  and  population  cumulative  distributions,  then  the  required 
sample  size  is  given  by  the  equation 

c 
<k(n)  -  /n"  "  &  <46> 

where  ca  is  the  constant  defining  the  critical  value  of  the  d  statistic 
in  the  Kolmogorov-Smirnov  test  for  samples  of  size  greater  than  35. 
Knowing  the  values  of  a  and  3,  it  is  then  possible  to  pick  out  the 
value  of  Cq^  from  statistical  tables  (Massey,  1951;  Ostle,  1964)  and 
hence  the  value  of  n,  being  the  requisite  number  of  samples  to  meet 
the  stipulated  criteria  could  be  computed. 

SIMULATION  IN  WATER  RESOURCES  SYSTEMS 

Simulation  means  to  duplicate  the  essence  of  the  system  or 
activity  without  actually  attaining  reality  itself.  Simulation  has 
been  used  traditionally  in  engineering.  The  use  of  conceptual  system 
models,  scale  models,  analogs,  and  laboratory  experimentation  are  but 
some  of  the  general  simulation  techniques  used  in  engineering.  However, 
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this  study  deals  with  the  simulation  of  processes  involved  in  the 
prediction  of  a  stream's  waste  assimilative  capacities  using  a 
digital  computer*  The  approximation,  complexity,  accuracy  and  cost  of 
simulation  are  somewhere  between  those  of  mathematical  analysis  and 
prototype  testing.  Simulation  is  especially  adaptable  to  the  study 
of  complex  systems  and  in  a  way,  it  may  be  considered  as  a  numerical 
method  for  the  solution  of  complicated  probabilistic  or  stochastic 
processes*  It  is  used  to  circumvent  the  difficulties  of  duplication 
of  environment,  of  mathematical  formulation,  of  lack  of  analytic 
solution  techniques,  or  of  experimental  impossibilities* 

Analog  and  digital  computers  are  frequently  used  for  pre- 
dicting the  responses  of  certain  water  quality  parameters  in  surface 
water  bodies  like  rivers,  lakes,  estuaries,  etc*  With  the  advent  of 
modern  electronic  computers,  more  and  more  efforts  are  being  directed 
toward  water  resources  planning  on  a  regional  basis,  with  optimization 
of  water  and  waste  water  treatment  processes  taken  as  a  whole  instead 
of  treating  these  two  as  separate  entities  as  has  been  done  in  the  past. 
Since  water  resources  systems  are  asually  complex  with  multiple  alter- 
natives for  design,  allocation,  and  operation,  simulation  seems  to  be 
a  very  useful  tool  in  water  resources  system  design* 

When  the  mathematical  model  for  a  process  has  been  decided 
upon,  various  elements  of  the  process  can  be  represented  in  the  computer 
so  that  the  outputs  of  one  part  of  the  system  constitute  the  inputs 
into  one  or  more  other  elements*  A  complex  system  may  thus  be  represented 
by  a  series  of  relatively  simple  elements.  In  such  a  system,  it  is 
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possible  to  vary  the  parameters  of  the  system  to  Investigate  some 
specific  aspects  of  the  system,  and  when  necessary,  to  adopt  systematic 
or  trial  and  error  procedures  for  the  optimization  of  the  system. 
Because  of  these  facilities,  digital  and  analog  computers  have  been 
used  extensively  for  the  simulation  of  dynamic  physical  systems. 

This  study  deals  with  the  simulation  of  a  river  system  using 
a  digital  computer  with  specific  reference  to  the  simulation  of  the 
phenomenon  of  a  stream's  assimilative  capacity  for  the  organic  waste 
discharged  into  it  without  creating  nuisance  conditions  in  the  stream. 
This  differs  from  analog  models  in  that  the  generated  data  are  used  for 
the  simulation  of  the  process  facilitating  the  study  of  the  probabil- 
istic aspects  of  the  process  as  a  whole.  Generated  data,  as  used  here, 
will  refer  to  generating  random  values  for  variables  in  the  process 
under  consideration  on  the  basis  of  statistical  models  for  the  vari- 
ables, which  could  not  be  distinguished  from  the  real  or  historical 
data  by  means  of  the  usual  statistical  tests  of  significance.  The  gen- 
erated data  did  not  occur  in  the  past,  nor  will  they  occur  in  the 
future;  but  based  on  certain  statistical  considerations,  it  could  have 
been  the  real  record.  Regarding  validity  of  simulation  Flagle  et  al. 
(1960)  state: 

•••It  must  be  remembered  that  simulation  yields  only 
an  empirical  form  of  knowledge,  fraught  with  the 
danger  that  the  stochastic  model  is  not  truly  repre- 
sentative. In  this  sense  it  is  admittedly  inferior 
to  mathematical  techniques  which  yield  functional 
relationship  between  the  variables.  Ultimately  we  hope 
that  the  more  slowly  developing  mathematical  capa- 
bilities will  make  simulation  as  we  practice  it,  un- 
necessary. Until  that  time  comes,  however,  simulation 
techniques  stand  as  powerful  aides  to  the  operations 
analyst  who  must  produce  useful  results. 
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MODEL  FOR  STREAM  ASSIMILATIVE  CAPACITY 
Factors  Affecting  DO  and  BOD 

Since  Streeter  and  Phelps  (1925)  propounded  their  theory 
defining  two  mechanisms  which  affect  the  DO  and  BOD  relationship 
namely,  the  BOD  removal  due  to  bacterial  oxidation  and  oxygen  addition 
due  to  atmospheric  reaeration,  several  other  factors  affecting  DO-BOD 
relationship  have  been  postulated  (Camp,  1963;  Dobbins,  1964;  Gannon, 
1963;  Thomas,  1948),  Some  of  the  major  factors  which  affect  the  as- 
similative capacity  of  a  stream  in  general  can  be  summarized  as  follows: 

1,  Removal  of  BOD  by  bacterial  oxidation  of  both  carbon- 
aceous and  nitrogenous  matter* 

2.  Removal  of  BOD  by  sedimentation  or  adsorption. 

3*  Addition  of  BOD  along  a  stretch  by  the  scour  of  bottom 
deposits  or  by  the  diffusion  of  partly  decomposed  organic  matter  from  the 
benthal  layer  into  the  overlying  water. 

4.  Addition  of  BOD  along  the  stretch  by  local  run  off. 

5.  Biological  removal  and  accumulation  of  BOD  by  fixed 
plants  and  algae. 

6.  Addition  of  oxygen  from  the  atmosphere. 

7.  Addition  of  oxygen  by  the  photosynthetic  action  of 
plankton  and  fixed  plants. 

8.  Removal  of  oxygen  by  the  respiration  of  plankton  and 
fixed  plants. 

Thomas  (1948)  introduced  a  rate  constant  K3,  with  units  of  per  day,  as 
a  means  of  accounting  for  the  removal  or  addition  of  BOD  by  deposition 
or  resuspension.  The  net  rate  was  assumed  to  be  proportional  to  K3L,  a 
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negative  value  of  K~  indicating  deposition  and  a  positive  value 
indicating  resus pens ion. 

In  considering  the  stabilization  of  organic  matter  in 
natural  streams,  distinction  should  be  made  between  those  factors 
which  affect  its  removal  without  necessarily  utilising  oxygen  and 
those  which  simultaneously  utilize  oxygen  while  removing  the  organic 
matter*  In  the  first  category  are  such  phenomena  as  sedimentation, 
scour*  and  biological  accumulation  and  the  second  category  includes 
various  forms  of  chemical  and  biological  oxidation.  Generally  when  a 
waste  is  discharged  to  a  stream*  it  is  extremely  difficult  to  isolate 
the  effects  due  to  each  of  these  factors  individually*  O'Connor  (1967) 
introduced  the  parameter  Kr*  with  units  of  per  day,  to  define  the  total 
rate  of  removal  of  BOD  in  the  river*  It  is  determined  from  a  series 
of  measurements  of  the  BOD  at  a  number  of  stations  downstream  from  the 
outfall*  Since  the  rate  of  removal  of  organic  matter  is  not  necessarily 
equal  to  the  rate  at  which  oxygen  is  utilized,  the  difference  between 
the  BOD  removal  rate  and  oxygen  utilization,  if  such  a  difference 
exists,  may  be  attributed  to  sedimentation,  scour,  benthal  demand,  and 
biological  accumulation*  If  none  of  these  factors  is  present  or  sig- 
nificant, it  is  usually  assumed  that  the  rate  at  which  the  organic 
matter  is  removed  is  equal  to  the  rate  at  which  the  dissolved  oxygen 
is  utilized  (O'Connor,  1967),  in  which  case  Kx  -  Kr. 

The  sources  of  oxygen  in  a  river  stretch  are  the  amount  of 
DO  in  the  Incoming  flow,  that  due  to  natural  or  artificial  reaeration 
and  that  due  to  the  photosynthetic  addition.  The  rate  of  reaeration 
is  proportional  to  the  dissolved  oxygen  deficit.  The  photosynthetic 
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source  depends  upon  many  factors  such  as  sunlight,  temperature,  mass 
of  algae  and  nutrients.  If  the  effect  of  these  factors  is  included  in 
the  term  Pfct,  representing  the  overall  rate  at  which  oxygen  is  released 
by  photosynthesis  and  if  the  photosynthetic  rate  if  assumed  to  vary  as 
the  sunlight,  reaching  a  peak  at  noon  and  a  zero  value  at  sunrise  and 
sunset,  then  this  source  could  be  defined  by  a  periodic  function 
(O'Connor,  1967)  as: 


P„  sin  St»   when  0  <  t*  <  p 

V  -   {  D    ^         -   -  (*7) 

0  when  1  -  p  <  t'  <  1 


where  p  is  the  period  of  oscillation  of  the  periodic  function,  Pm 
is  the  amplitude  of  the  periodic  function  and  t'  is  the  time  of  flow 
in  days  of  a  mass  of  water  since  sunrise  after  the  mass  of  water  has 
entered  the  head  end  of  a  river  reach  considered.  If  the  period  p  is 
twelve  hours  or  half  a  day,  the  periodic  function  shown  in  Eq,  47 
reduces  to 


P  sin  2nt*   when  0  <  t«  <  1/2 
Pt,  -   {  (47a) 

0  when  1/2  <  t»  <  1 


This  periodic  function  can  be  described  by  a  Fourier  series  as 
defined  below: 


00  00 

Pt«  -  a0  +   T      a^  cos  nt1  +   £   bn  sin  nt*       (48) 
n  -  1  n  -  1 
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The  coefficients  of  the  Fourier  series  can  be  evaluated  using 
Euler's  formula  (Kreyzig,  1962)  giving 


Pm 

*o  -  T  <49> 


an  "  fn^^tS^^  "W  4  I  £*  t   *>      <50) 


When  n  is  odd 


a„  -  0  (50a) 

n 


and  when  n  is  even 

2P 


m 


Also 


(n  -  l)(n  +  1)tt 


pm 


(50b) 


b!  -  Y  (51) 


and 


b_  -  0  for  n  -  2,  3,...  (51a) 


Hence 


P    P  2P 

Pt,  -  —  +  Y  sin   2TTt«  -  —  (I  cos  4TTt»  +  Ij  Snt*  +  ...)   (52) 

Thus  it  is  seen  that  the  oxygen  addition  due  to  photosynthesis  can  be 
represented  by  the  first  three  terms  of  Eq.  52  with  good  approximation. 
This  approach  of  expressing  quantitatively  the  diurnal  variations  in 
phot  ©synthetic  activity  is  decidedly  superior  to  averaging  over  the 
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entire  period  including  the  nighttime  flow  conditions  as  was  done 
by  Camp  (1963)  and  Dobbins  (1964), 

Bauat ions  for  BOD  and  DO  Profiles 

The  equations  for  the  BOD  and  DO  profiles  along  a  river 
stretch  are  based  on  the  following  assumptions: 

1.  The  stream  flow  is  steady  and  uniform, 

2.  Removal  of  BOD  by  both  bacterial  oxidation,  sedimentation 
or  adsorption  or  both  and  biological  accumulation  are  first  order  re- 
actions, the  rates  of  removal  being  proportional  to  the  amount  present. 

3.  The  removal  of  oxygen  by  benthal  demand,  by  plant 
respiration  and  the  addition  of  BOD  from  the  benthal  layer  are  uniform 
along  the  stretch. 

4.  Diurnal  variations  of  oxygen  addition  due  to  photosyn- 
thesis can  be  represented  by  a  sine  function. 

5.  The  rate  coefficients  affecting  the  oxidation  of  organic 
matter  and  atmospheric  reaeration  can  be  taken  as  random  variables  and 
these  could  be  defined  by  suitable  statistical  models. 

6.  The  BOD  and  DO  are  uniformly  distributed  over  each  cross 
section  so  that  the  equations  can  be  written  in  the  one-dimensional 
form. 

Under  the  foregoing  assumptions,  the  differential  equation 
for  the  BOD  profile  is  given  by 


f  -  -  «i  +  K3  *  VL  (53) 
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where  Kfe  is  the  coefficient  for  determining  the  rate  of  BOD  removal 
due  to  biological  accumulation.  The  differential  equation  for  the 
oxygen  profile  is  given  by 

~  -  Kxh   -  K2D  -  Pffl  (~  +  |  sin  2-nt'    -  ^  cos  4TTt')  +  R  +  Db   (54) 

where  R  is  the  rate  of  oxygen  removal  by  algal  respiration  and  IX  It  the 
benthal  demand  both  with  the  units  of  mg/(l)(day).  When  these  two 
equations,  Eqs.  53  and  54,  are  solved  simultaneously  with  known  initial 
conditions  for  the  BOD  and  DO  at  the  upper  end  of  a  river  section,  it 
is  possible  to  predict  the  state  of  DO  at  any  point  within  the  reach. 

Application  to  the  Problem  under  Study 

Generated  data  for  the  values  of  K.,  K2  and  DO  found  necessary 
for  the  probabilistic  analysis  of  the  process  under  investigation  were 
developed  by  using  Monte  Carlo  methods  and  simulation  techniques.  The 
statistical  models  as  developed  in  this  study  for  K.  and  K2  were  used 
in  generating  data  applicable  to  the  Cincinnati  Pool  of  the  Ohio  River 
prior  to  the  construction  of  Markland  Dam.  Parameters  pertaining  to 
other  mechanisms  namely  photosynthesis,  respiration,  etc.,  are  evalu- 
ated from  the  published  Ohio  River  Survey  data  (USPHS,  1960).  These 
parameters  were  taken  at  their  average  values  and  the  expected  and 
medal  values  for  DO  as  predicted  by  the  probabilistic  model  for  known 
initial  conditions  of  BOD  and  DO  deficit  were  compared  with  the  actual 
DO  concentrations  observed  in  the  survey. 


94 


Using  a  subroutine  available  at  the  Digital  Computer  Labora- 
tory at  the  University  of  Illinois,  two  sets  of  uniformly  distributed 
random  number  sequences  were  generated  and  transformed  to  correspond 
to  normally  distributed  random  variables  with  given  characteristics. 
The  first  set  of  generated  values  was  used  to  evaluate  K^  and  the 
second  K2,   In  the  case  of  K1#  the  generated  normally  distributed 
values  with  parameters  matching  the  historical  values  themselves  con- 
stitute the  simulated  values  for  K, .  In  the  case  of  K2,  the  generated 
values  constitute  the  random  components  of  K->  and  when  these  are  added 
to  the  trend  component  as  given  by  the  regression  equation,  Eq.  31, 
the  random  variable  K2  is  obtained*  The  generated  data  were  tested  to 
see  whether  the  sample  could  be  distinguished  from  the  theoretical 
model  using  the  Kolmogorov-Smirnov  test.  The  generated  values  for  K^ 
and  K2  along  with  other  parameters  as  determined  in  the  Survey  (USPHS, 
1960)  were  used  in  the  simulation  model  to  obtain  probabilistic  estimates 
of  DO  at  the  desired  sections  in  the  river  reach  considered. 

It  is  appropriate  to  quote  what  Flagle  et  al.  (1960)  had  to 

say  on  probabilistic  approaches  for  solving  problems  which  do  not 

strictly  conform  to  deterministic  rules: 

...as  soon  as  human  participation  in  the  system  operation 
occurs  it  is  necessary  to  introduce  arguments  of  proba- 
bility and  calculate,  not  how  the  system  will  certainly 
and  invariably  perform,  but  only  how  it  will  perform  a 
certain  fraction  of  the  time  it  is  stated.  Sometimes, 
as  in  the  case  of  the  atomic  bomb  already  mentioned, 
statistical  methods  must  be  used  even  for  a  mechanical 
system  when  the  laws  of  physics  that  apply  (like  those 
for  radio-activity)  are  inherently  statistical  in  nature. 

This  study  is  limited  to  the  development  and  subsequent 

testing  of  the  generated  K^  and  K2  data  and  the  analysis  of  some  of 
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the  probabilistic  characteristics  of  the  DO  response  in  stream  as- 
similative capacity  determinations*  Future  studies  may  deal  with 
the  simulation  of  systems  of  waste  treatment  facilities,  and  the 
combination  of  water  and  waste  water  treatment  facilities  for 
determination  of  the  optimal  design. 
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VI.   APPLICATION  OF  MONTE  CARLO  METHOD 
TO  THE  OHIO  RIVER  SURVEY  DATA 


CHOICE  OF  PUBLISHED  RIVER  SURVEY  DATA 

The  primary  objective  of  this  study  is  to  develop  a  method 
for  estimating  the  stream's  waste  assimilative  capacity  taking  into 
account  the  variability  of  the  reaction  coefficients  affecting  the 
DO-BOD  relationship;  the  need  for  which  has  long  been  felt.   In  order 
to  express  quantitatively  the  variability  in  these  coefficients 
using  probability  measure,  extensive  data  are  needed  and  published 
data  for  five  river  surveys  (Gannon,  1963;  The  Resources  Agency  of 
California,  1962;  TVA,  1962;  USPHS,  1960;  USPHS,  1963)  were  examined 
to  study  in  detail  the  probabilistic  variations  in  these  coefficients 
and  to  verify  the  observed  data  of  dissolved  oxygen  with  the  predicted 
values  using  Monte  Carlo  methods.  Though  the  Sacramento  River  survey 
(The  Resources  Agency  of  California,  1962)  was  extensive,  encompassing 
information  on  physical,  chemical,  and  biological  aspects,  photosyn- 
thesis, respiration,  etc.,  the  number  of  observations  made  for  long- 
term  BOD  progression  below  a  single  major  waste  discharge  was  far  less 
than  that  made  in  the  Ohio  River  survey  (USPHS,  1960),  which  was  other- 
wise comparable  to  Sacramento  River  survey  in  various  other  aspects 
observed.  The  long-term  BOD  experiment  results  published  by  Gannon 
(1963)  for  the  Clinton  and  Tittabawassee  Rivers  were  carried  out  at 
different  conditions  of  incubation  temperature,  mixing,  dilution  and 
nitrification  inhibition  agents;  whereas  in  the  Ohio  River  survey  tests 
for  the  BOD  progression  river  samples  were  all  carried  out  according 
to  the  procedures  stipulated  in  Standard  Methods  (APHA,  1955).  Also 
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no  attempts  were  made  in  this  survey  (Gannon,  1963)  to  evaluate  the 
effect  of  photosynthesis  due  to  aquatic  vegetation,  even  though  the 
author  concluded  that  the  abundance  of  aquatic  vegetation  was  respon- 
sible for  a  manifold  increase  in  the  river  BOD  removal  rate  compared 
to  the  BOD  removal  rate  observed  in  bottle  experiments.  The  report 
on  the  Illinois  River  System  (USPHS,  1963)  did  not  include  details  of 
the  BOD  progression  data  and  the  number  of  experiments  carried  out  in 
this  study  to  evaluate  K,  values  was  too  few  to  make  reliable  statis- 
tical inference.  Hence  it  was  observed  that  the  data  published  for 
the  Ohio  River-Cincinnati  Pool  survey  (USPHS,  1960)  which  included 
long-term  BOD  results  on  river  samples,  dark  and  light  bottle  tests 
for  determining  photosynthesis  and  respiration  of  phy top lank tons, 
benthal  deposits,  hydraulic  characteristics  of  the  river,  etc*,  are 
best  suited  for  this  study. 

In  any  river  survey  for  determining  the  pollution  assimilative 
capacity,  it  is  necessary  to  determine  the  river  reaeration  coefficient, 
K2,  independently.  In  order  to  estimate  the  random  values  of  K2  for 
the  Ohio  River-Cincinnati  Pool  reach,  it  is  assumed  that  the  regression 
equation  for  the  trend  component  and  the  probability  distribution 
function  Cor  the  random  component  of  K2  as  developed  for  Tennessee 
Valley  rivers  are  applicable  to  the  Ohio  River  reach  under  consideration. 

VALUES  OF  THE  PARAMETERS  USED  IN  THE  MODEL 
Deoxygenation  Coefficient 

With  the  advent  and  increased  use  of  secondary  treatment 
facilities  for  waste  disposal,  the  importance  of  considering  oxygen 
demand  due  to  nitrification  has  been  recognized.  O'Connor  (1967), 


98 

and  Stratton  and  McCarty  (1967)  proposed  mathematical  models  to  account 
for  the  nitrogenous  oxygen  demand  in  oxygen  balance  studies  in  rivers. 
As  indicated  in  Chapter  IV,  in  the  case  of  the  Ohio  River  study,  neither 
the  treatment  plant  effluents  nor  the  river  samples  obtained  below  the 
waste  outfalls  showed  any  significant  nitrification  for  about  10  days. 
Though  the  mathematical  model  proposed  to  be  used  for  the  Ohio  River 
does  not  contain  terms  to  account  for  nitrification,  conceptually  it 
does  not  present  any  difficulty  to  include  terms  to  account  for  the 
oxygen  consumption  due  to  nitrification  for  situations  where  this  ef- 
fect is  significant.  Thus,  the  only  factor  which  affects  BOD  removal 
with  concomitant  removal  of  oxygen  is  that  due  to  the  biological  oxi- 
dation of  organic  matter.  The  details  of  the  analysis  of  BOD  progression 
data  obtained  on  Ohio  River-Cincinnati  Pool  samples  were  presented  in 
Chapter  IV,  The  variations  of  K^  values  for  these  river  samples  were 
found  to  be  random  and  adequately  defined  by  a  statistical  model)  namely, 
that  the  Kj^  values  for  river  samples  were  normally  distributed  with  a 
mean  of  0,173  per  day  and  a  standard  deviation  of  0.066  per  day  having 
a  coefficient  of  variation  of  38  percent.  This  statistical  model  for 
Kj  is  adapted  for  simulation  studies  using  the  Monte  Carlo  method. 

Effect  of  Sludge  Deposits 

One  of  the  greatest  difficulties  in  stream  surveys  for 
estimating  assimilative  capacities  is  to  evaluate  the  values  of  K-j, 
K5,  and  Dfe  individually.  To  overcome  the  difficulty  of  estimating  K3 
and  Kb  values  Velz  and  Gannon  (1962),  Gannon  (1963),  O'Connor  (1967), 
and  others  advocate  the  use  of  the  term  Kr  which  defines  the  coefficient 
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for  the  BOD  removal  rate  due  to  the  lumped  effects  of  K^,  K3,  K., 
and  other  unaccounted  causes.   In  the  absence  of  significant  sludge 
deposits  and  fixed  vegetation,  the  values  of  the  rate  coefficients  K3 
and  Kb  can  be  taken  as  zero  leaving  only  the  oxidation  of  organic 
matter  as  the  significant  operative  mechanism  responsible  for  the 
removal  of  BOD  from  the  system.  Bottom  samples  collected  during  the 
1957  study  of  the  Ohio  River-Cincinnati  Pool  (USPHS,  1960)  revealed 
that  there  rat  no  significant  sludge  deposits  in  the  reaches  in- 
vestigated to  warrant  recognition  of  this  as  an  important  factor  in 
the  oxygen  balance  studies.  This  is  further  substantiated  by  the 
experience  of  several  investigators  (Nejedly,  1966;  Thomas,  1948; 
Velz,  1958),  that  the  problem  of  sludge  deposits  is  greatly  allevi- 
ated if  settleable  solids  are  removed  from  the  raw  waste  prior  to 
discharge  into  rivers.  Since  all  the  wastes  emanating  from  Cincinnati, 
Ohio  received  primary  treatment,  and  since  the  survey  report  indicates 
that  the  sludge  deposits  were  insignificant,  all  the  mechanisms  of  BOD 
and  DO  removal  governed  by  sludge  deposits  can  be  treated  as  insig- 
nificant and  consequently  K~  can  be  taken  as  zero. 

In  deep  rivers  such  as  the  Ohio  River,  BOD  removal  due  to 
extraction  and  accumulation  of  fixed  plants  is  insignificant  (Gannon, 
1963;  Heukelekian,  1967).   It  is  relevant  here  to  quote  what  Gannon 
(1963)  has  to  say  on  river  BOD  removal  mechanisms: 

Not  all  rivers  have  high  BOD  removal  and  many  have 
rates  which  closely  parallel  the  rates  determined  in 
the  laboratory  bottle  experiments.  This  appears  to 
be  particularly  true  for  the  large  rivers,  where  there 
is  relatively  less  contact  with  river  bottom  and  sides 
and  where  there  is  no  gross  dispersed  type  of  growth. 
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Thus,  it  is  seen  that  the  only  significant  mechanism  of  any  consequence 
affecting  the  BOD  removal  in  the  Ohio  River-Cincinnati  Pool  area  is 
that  due  to  the  biological  oxidation  of  organic  matter.  Also  con- 
sequent to  the  absence  of  significant  sludge  deposit,  the  benthal 
demand  for  oxygen  is  negligible, 

Reaeration  Coefficient 

Churchill  et  al.  (1962)  suggest  that  when  the  regression 
equation,  Eq.  23,  is  applied  to  polluted  streams,  the  basic  reaeration 
rates  should  be  modified  on  a  percentage  basis  as  indicated  by  rela- 
tive reaeration  rates  determined  in  the  laboratory  by  experiments  for 
the  polluted  water  and  for  water  samples  obtained  upstream  of  waste 
discharge.  It  is  generally  considered  that  the  constituents  in 
waste  waters  tend  to  reduce  reaeration  rates.  Downing  et  al.  (1957), 
and  Downing  and  Truesdale  (1955)  conducted  experiments  to  determine 
the  effects  of  several  contaminants  and  mixtures  of  contaminants  on 
the  exchange  coefficient  for  oxygen  in  water  agitated  at  different 
rates  in  laboratory  absorption  vessels.  They  concluded  that  house- 
hold detergents  reduced  exchange  coefficients  in  clean  water  by 
amounts  which  depended  on  their  initial  concentration  and  rate  of 
agitation.  They  further  concluded  that  the  effects  of  settled  sewage 
in  the  absence  of  any  added  anionic  detergent  were  lower  than  when 
these  materials  were  present.  The  reported  percent  reduction  in  K2 
values,  under  laboratory  experimental  conditions  due  to  municipal 
wastes  varies  from  10  to  30  depending  on  the  type  and  concentration 
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of  the  constituents  and  the  rate  of  agitation  (O'Connor,  1958; 
Downing  et  al,,  1957;  Poon  and  Campbell,  1967)*  Poon  and  Campbell 
(1967)  found  that  at  low  concentrations  of  suspended  solids  in  tap 
water,  the  transfer  rates  were  enhanced  to  the  extent  of  10  to  30 
percent*  Since  the  characteristics  and  concentrations  of  waste 
discharges  are  subject  to  considerable  fluctuations,  no  attempt  has 
been  made  to  assign  a  numerical  value  for  the  possible  effect  of 
the  pollutants  on  the  mean  river  reaeration  coefficient  as  predicted 
by  the  regression  equation  used  in  this  study.  Since  the  reaeration 
coefficient  is  treated  as  a  random  variable  comprising  of  trend  and 
random  components,  the  error  due  to  neglecting  the  effects  of  waste 
constituents  which  have  compensating  effects  on  the  reaeration  coef- 
ficient is  considered  negligible. 

The  values  of  velocity  of  flow,  and  depth  used  for  deter- 
mining the  trend  components  of  K2  are  shown  in  Table  11,  for  six 
different  cases  in  which  observations  were  made  at  two  different 
sections  of  the  river  reach  under  consideration.   In  these  cases, 
the  downstream  samples  were  obtained  presumably  after  a  time  lag  equal 
to  the  flow-through  time  from  the  upstream  section  of  the  reach  with 
the  result  that  the  same  body  of  water  had  been  sampled  as  it  flowed 
down  the  river.  The  percent  deviation  of  actual  K2  value  from  the 
predicted  value  was  taken  to  be  normally  distributed  with  mean  zero 
and  standard  deviation  of  36,8  percent.  The  order  of  magnitude  of  the 
residual  error  distribution  of  K~  values  applicable  to  the  Ohio  River 
study  has  been  assumed  to  be  the  same  as  that  for  the  Tennessee  Valley 
rivers. 
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Photosynthesis  and  Respiration 

Values  for  these  parameters  computed  from  the  dark  and  light 
bottle  observations  reported  for  the  Ohio  River  (USPHS,  1960)  are 
shown  in  Table  11,  Since  photosynthesis  and  respiration  of  phyto- 
plankton  are  greatly  dependent  on  teaperature  and  sunlight  intensity, 
it  is  likely  that  these  parameters,  over  a  long  period  of  observation 
will  tend  to  show  a  trend  in  their  values.  It  is  more  realistic  to 
consider  these  parameters  as  variables  with  trend  and  random  components 
since  the  factors  which  affect  these  parameters  fluctuate  having  a 
trend  with  diurnal  variations,  Siace  the  scope  of  the  Ohio  River  sur- 
vey light  and  dark  bottle  studies  was  limited  to  12  days  spread  over 
two  months,  there  is  not  enough  information  to  discern  these  trends. 
As  the  maximum  values  for  phot ©synthetic  oxygen  addition  are  assumed 
to  occur  around  mid-day,  the  number  of  such  observations  made  during 
the  survey  Is   inadequate  to  formulate  and  test  any  statistical  model 
to  characterize  these  parameters.  Only  arithmetic  average  values  of 
these  parameters  obtained  from  the  observations  for  different  periods 
are  used  in  this  study.  It  is  conceptually  feasible  to  consider  these 
two  parameters  also  as  random  variables  provided  enough  data  could  be 
collected  to  establish  the  trend  and  random  components. 

SAMPLE  SIZE  IN  SIMULATION  STUDIES  OF  THE  OHIO  RIVER 

It  was  indicated  earlier  that  the  required  sample  size  in 
simulation  studies  could  be  determined  using  Eq.  46,  if  the  desired 
precision  and  confidence  levels  are  known.  Flagle  (1960)  suggests 
values  of  10  percent  and  95  percent  for  error  and  confidence  levels 
respectively.  Referring  to  Table  1  of  Massey  (1951)  for  the  critical 
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values  of  c^(n)  of  the  maximum  absolute  difference  between  sample  and 
population  cumulative  distributions,  the  critical  value  of  'd'  statistic 
in  the  Kolmogorov-Smirnov  test  for  95  percent  confidence  level  is  given 


"0.05  "  '-f- 
/n 


Thus,  substituting  the  values  for  c^  and  3,  being  1,36  and  10  percent 
respectively  in  Eq.  46,  and  solving  for  n,  the  required  sample  size 
is  indicated  as  185  or  approximately  200. 

Hence  In  order  to  have  confidence  level  of  95  percent  that 
the  information  on  required  distribution  is  not  different  from  the  actual 
value  by  more  than  10  percent,  the  sample  size  is  to  be  200 . 

RESULTS  OF  THE  MONTE  CARLO  METHOD  APPLIED  TO  THE  OHIO  RIVER 

A  flow  diagram  for  the  computer  simulation  studies  using  the 
Monte  Carlo  method  following  the  general  procedure  enumerated  in  Chapter 
V  is  shown  in  Appendix  C.  Two  separate  sets,  each  having  200  random 
numbers  with  uniform  distribution  (0,1)  were  generaged  using  a  sub- 
routine available  with  the  University  of  Illinois  Digital  Computer 
Laboratory,  The  first  set  of  random  numbers  was  transformed  to  corres- 
pond to  K.  values  with  a  mean  value  of  0.173  per  day  and  standard 
deviation  of  0.066  per  day.  The  second  set  of  random  numbers  was 
transformed  to  correspond  to  the  random  variations  in  percent  error  of 
predicted  K2  values  with  a  mean  value  of  zero  and  standard  deviation 
of  38,6  percent. 
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It  is  likely  that  some  extremely  and  unreasonably  high  or 
low  values  will  be  generated  by  the  Monte  Carlo  method  and  it  becomes 
necessary  to  ignore  such  values  (Montgomery  and  Lynn,  1964)  or  these 
have  to  be  corrected  on  a  reasonable  basis  (Ramaseshan,  1964),  In 
this  study  the  transformed  variables  with  values  beyond  twice  the 
deviation  from  the  stipulated  mean  were  assigned  mean  values.  The 
generated  random  variates  were  tested  for  goodness  of  fit  with  the 
assumed  distributions  using  the  Kolmogorov-Smirnov  test.  The  distri- 
butions of  the  generated  values  were  found  to  satisfy  the  test. 

With  known  conditions  for  velocity  and  depth  of  flow,  the 
trend  component  of  K2  could  be  evaluated  using  Eq.  31,  Knowing  the 
mean  K2  values  for  the  given  conditions  of  flow,  the  random  variations 
in  K2  could  be  computed  from  the  generated  random  variations  in  percent 
error.  These  generated  random  values  for  the  variations  in  K2  values 
were  added  algebraically  to  the  mean  value  to  obtain  the  random  values 
of  K2, 

The  river  system  for  determining  the  assimilative  capacity 
was  simulated  using  Eqs,  53  and  54  in  which  K3,  K,,  and  D.  were  taken 
as  zero  for  reasons  discussed  earlier.  Generated  values  of  Kj,  K2, 
along  with  mean  values  for  Pm  and  R  indicated  in  Table  11  were  used  in 
determining  the  response  of  the  system.  The  values  for  DO  deficit 
predicted  by  Monte  Carlo  techniques  were  verified  with  observed  values 
in  six  different  cases  in  which  the  same  body  of  river  water  had  been 
sampled  at  two  different  sections  downstream  of  all  the  major  waste 
discharges  emanating  from  Cincinnati,  Ohio,  In  each  of  these  cases, 

the  process  was  simulated  with  known  initial  conditions  and  the  state 

i 

j 
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of  dissolved  oxygen  in  terms  of  oxygen  deficit  was  predicted.  The 
initial  conditions  of  ultimate  BOD,  dissolved  oxygen,  temperature, 
and  the  final  conditions  of  dissolved  oxygen  and  tejiperature  along 
with  the  flow-through  time  between  the  two  sections  are  presented  in 
Table  12,  for  each  of  the  six  cases  considered. 

Two  hundred  values  for  dissolved  oxygen  deficit  values  were 
generated  for  each  case  study  solving  the  differential  equations  using 
a  subroutine  in  the  University  of  Illinois  Digital  Computer  Laboratory. 
The  subroutine  employs  fourth  order  Runge-Kutta  method  for  solving  dif- 
ferential equations,  the  details  of  which  can  be  found  in  standard 
textbooks  on  numerical  analysis  (Fox,  1962;  McCracken  and  Dorn,  1964). 
Utilizing  the  generated  values  of  DO  deficit,  the  frequency  distribution 
of  probable  values  was  computed  and  these  are  shown  in  Table  13,  and 
plotted  in  Figures  14  through  19.  The  expected  and  most  probable  values 
of  DO  deficits  for  each  of  these  cases  are  shown  in  Table  14.  In  order 
to  evaluate  and  compare  the  results  obtained  by  Monte  Carlo  techniques, 
predicted  values  of  DO  deficits  were  computed  for  each  case  using  equa- 
tions postulated  by  Streeter  and  Phelps  (1925),  and  Camp  (1963).  Average 
value  of  0.173  per  day  for  K^  and  the  value  obtained  from  the  regression 
equation,  Eq,  31,  for  K2  using  the  hydraulic  characteristics  reported 
in  Table  11  were  used  in  these  formulae.  These  results  are  presented 
in  Table  14,  Percent  error  in  predicting  DO  deficit  was  computed  by 
multiplying  100  with  the  fraction  obtained  by  dividing  the  absolute 
difference  between  actual  and  predicted  values  of  DO  deficit  with  the 
actual  deficit. 
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FIGURE  14.     FREQUENCY  DISTRIBUTION  OF  DO  DEFICIT 
IN  CASE  STUDY  NO.   1 
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FIGURE  15.     FREQUENCY  DISTRIBUTION  OF  DO  DEFICIT 
IN  CASE  STUDY  NO.  2 
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FIGURE  16.     FREQUENCY  DISTRIBUTION  OF  DO  DEFICIT 
IN  CASE  STUDY  NO.   3 
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FIGURE  17.     FREQUENCY  DISTRIBUTION  OF  DO  DEFICIT 
IN  CASE  STUDY  NO.  4 
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FIGURE  18.  FREQUENCY  DISTRIBUTION  OF  DO  DEFICIT 
IN  CASE  STUDY  NO.  5 
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FIGURE  19.     FREQUENCY  DISTRIBUTION  OF  DO  DEFICIT 
IN  CASE  STUDY  NO.  6 
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VII.   DISCUSSION 

The  significance  and  importance  of  considering  the  variations 
in  K  and  K  were  brought  out  in  Chapter  III.  In  the  carefully  selected 
river  reaches  of  the  Tennessee  Valley  rivers,  the  atmospheric  reaeration 
was  found  to  vary  significantly  in  every  reach  considered,  even  though 
each  set  of  observations  was  made  under  constant  river  flow  conditions. 
The  observed  range  of  values  for  K2  was  considerably  more  for  shallow 
and  high  velocity  flows  than  for  deep  and  low  velocity  flows.  As 
pointed  out  earlier,  in  one  set  of  experiments  in  the  Holston  River, 
K2  had  a  range  of  values  from  0.10  to  1.18  per  day  with  an  arithmetic 
average  of  0.63  per  day.  All  the  mathematical  models  developed  for 
predicting  the  reaeration  rate  coefficient  attempt  to  estimate  only 
the  average  value  for  K2.  Efforts  so  far  made  to  evaluate  the  stream 
assimilative  capacities  fail  to  take  into  account  the  variations  in 
the  actual  values  of  K2  from  the  predicted  values.  These  variations 
are  considerable  and  lead  to  significant  errors  in  predicting  the  DO 
response,  if  omitted. 

Attempts  have  been  made  to  relate  the  reaeration  coefficient 
to  such  hydraulic  characteristics  as  velocity  and  depth  of  flow,  energy 
slope,  surface  renewal  rates,  etc.  These  were  not  found  to  explain 
adequately  the  variations  in  the  predicted  values  for  K~.   Inclusion  of 
the  variables,  such  as  molecular  diffusion,  dispersion  coefficient, 
surface  tension,  etc.,  did  not  significantly  improve  the  predictions 
for  reaeration  coefficients  (TVA,  1962).  Variations  in  the  values  for 
reaeration  coefficient  may  be  attributed  to  the  errors  in  estimation 
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of  the  parameters  involved  in  the  predictor  equations,  simplif i cat ion 
in  the  assumptions  for  developing  the  equations  and  other  uncontrollable 
and  unaccounted  causes  such  as  wind  effects,  turbulence,  etc.  As  these 
factors  are  not  quantifiable  in  a  deterministic  sense,  the  combined 
effects  have  to  be  expressed  in  terms  of  probability  measure.  For  the 
Tennessee  Valley  rivers,  the  variations  in  reaeration  coefficients  have 
been  found  to  be  adequately  defined  by  Gaussian  distribution. 

Though  the  deoxygenation  and  reaeration  coefficients  were 
the  only  two  parameters  treated  as  random  variables  in  this  study, 
the  probabilistic  concepts  should  be  extended  to  other  factors  such  as 
photosynthesis,  respiration,  biological  extraction  of  BOD,  etc.,  which 
are  subject  to  random  fluctuations.  Also  the  covariance  of  these  factors 
needs  to  be  given  consideration.  Attempts  have  been  initiated  at  the 
University  of  Illinois  to  study  the  possible  effects  of  different  rates 
in  deoxygenation  on  atmospheric  reaeration  capabilities  in  laboratory 
channels.  Considerable  additional  research  is  needed  in  this  area. 

Monte  Carlo  techniques  applied  to  the  problem  of  determining 
the  DO  response  of  streams  receiving  organic  waste  loads  appear  to 
yield  satisfactory  results.  The  results  of  the  analysis  of  the  Ohio 
River  survey  data  presented  in  Table  14  indicate  that  there  is  con- 
siderable difference  between  the  DO  deficits  actually  observed  and  those 
predicted  by  the  deterministic  equations.  Application  of  the  Streeter- 
Phelps  equation  in  which  the  effects  of  photosynthesis  and  respiration 
are  omitted  and  in  which  the  reaction  coefficients  are  taken  as  constants 
in  a  river  reach,  results  in  prediction  errors  ranging  from  0.5  to  33.5 
percent.  Use  of  Camp's  equation  (1963)  which  accounts  for  an  average 
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value  of  photosynthesis  also  results  in  prediction  errors  ranging  from 
6.5  to  34,7  perdent.  The  use  of  an  average  value  for  photosynthesis 
in  Camp's  model  irrespective  of  daytime  or  nighttime  flow  conditions 
is  unrealistic  since  no  photosynthesis  can  be  expected  to  occur  during 
nighttime.  As  the  characteristics  of  the  waste  emanating  even  from  a 
single  source  vary  from  time  to  time  and  the  type  and  distribution  of 
microflora  stabilising  the  waste  organic  matter  are  bound  to  change 
in  the  natural  environment,  it  is  realistic  to  consider  all  the  possible 
variations  in  parameters  characterizing  these  changes  instead  of  as- 
suming them  to  be  time  invariant. 

The  percent  error  involved  in  predicting  DO  deficits  using 
the  probabilistic  model  range  from  0.4  to  15.7  when  the  predictions 
are  based  on  expected  values  and  range  from  0.0  to  11.5  when  based  on 
■est  probable  values.  It  is  seen  that  under  different  conditions,  the 
values  of  DO  deficit  predicted  by  the  probabilistic  model,  particularly 
the  most  probable  values,  tend  to  be  closer  to  the  observed  values  all 
the  time  than  the  values  predicted  by  the  deterministic  equations. 
The  ooncept  of  most  probable  value  is  not  unfamiliar  in  sanitary  engi- 
neering practice,  since  the  enumeration  of  coliform  organisms  in  water 
samples  involves  this  concept. 

The  average  value  for  the  percent  errors  in  DO  deficits  pre- 
dicted by  the  Streeter-Phelps  equation  for  the  six  cases  investigated 
in  this  study  is  14.0  and  the  corresponding  value  for  the  Camp  equation 
is  15.2.  The  average  values  in  percent  errors  using  the  probabilistic 
model  are  8,7  and  8.2  respectively  for  the  expected  value  and  most  prob- 
able value  cases.  Thus  it  is  seen  that  the  use  of  most  probable  value 
predicted  by  the  probabilistic  model  is  likely  to  give  more  reliable 
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information  on  the  DO  response  than  the  deterministic  equations 
considered  in  this  study.  However  the  Streeter-Phelps  model  has  pre- 
dicted the  DO  deficit  within  one  percent  error  when  the  initial  and 
final  observations  were  carried  out  at  night  (Case  Study  No.  2)  and 
the  error  in  prediction  using  Camp's  equation  for  this  case  is  6.5 
percent  which  is  the  best  value  obtained  employing  Camp's  model. 

The  errors  inherent  in  the  predictions  using  probabilistic 
model  should  be  considered  as  being  within  practical  limits,  since 
several  approximations  and  uncertainties  creep  in  the  evaluation  of 
various  factors  which  affect  the  DO-BOD  relationship  in  streams.  The 
estimation  of  time  of  flow  between  two  sections  in  a  river  reach  could 
at  best  be  only  approximate,  so  also  the  average  values  of  depths  and 
velocities  of  flow.  Any  errors  in  the  assumption  of  oxygen  saturation 
values  at  different  temperatures  are  bound  to  be  reflected  in  the  final 
results  of  DO  deficit  predictions.  Again  there  is  the  raging  controversy 
about  the  adequacy  of  BOD  bottle  experiments  to  duplicate  the  oxidation 
process  in  the  unrestricted  natural  environment.  Sanitary  engineers 
have  to  put  up  with  this,  being  the  only  expedient  method  available  at 
present  to  evaluate  the  biodegradability  of  organic  matter. 

The  use  of  the  regression  equation,  Eq.  31,  developed  from 
the  data  collected  in  Tennessee  Valley  rivers,  for  computing  mean 
values  of  K2  applicable  to  the  Ohio  River  involves  extrapolation  to 
conditions  other  than  those  for  which  the  equation  was  developed.  The 
validity  and  usefulness  of  the  equations  developed  by  O'Connor  and 
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Dobbins  (1956)  and  Dobbins  (1964)  were  not  established  beyond  reasonable 
doubt.  The  discrepancies  in  observed  and  predicted  values  using  these 
equations  have  been  pointed  out  (Churchill  et  al.,  1962;  Thackston  and 
Krenkel,  1965).  The  use  of  the  regression  equation  developed  for 
Tennessee  Valley  rivers,  being  the  best  available  tool  for  estimating 
the  reaeration  coefficient  independently,  might  be  another  source  of 
error,  since  this  involves  extrapolation  as  indicated  earlier. 

The  strength  of  the  proposed  method  lies  in  the  fact  that  in 
the  face  of  several  ambiguities,  it  predicts  a  range  of  values  for  DO 
deficits  with  associated  probabilities  instead  of  a  single  value  as  in 
deterministic  equations.  Thus  it  enables  one  to  quantify  the  uncer- 
tainties in  terms  of  probability  measure  and  to  consider  the  probability 
of  river  DO  being  equal  to  or  less  than  certain  concentrations  under 
given  waste  and  stream  flow  conditions. 

Figure  20  shows  the  probable  range  of  values  for  the  oxygen 
sag  profile  in  the  Ohio  River  and  the  associated  probabilities  at  a 
few  selected  sections  for  the  initial  conditions  represented  in  Case 
Study  No.  5.  The  histograms  are  drawn  with  class  intervals  of  0.25 
mg/1.  Class  intervals  smaller  than  0,25  mg/1  are  considered  unwarranted 
in  view  of  the  sensitivity  of  the  DO  determination  methods.  The  pro- 
files predicted  by  Camp's  and  Streeter-Phelps '  equation  are  also  shown 
in  Figure  20,  For  short  times  of  flow,  the  values  of  DO  predicted  by 
the  two  deterministic  equations  and  the  modal  value  of  the  probabilistic 
model  are  in  close  agreement  within  practical  limits.  As  the  time  of 
flow  increases,  there  is  considerable  deviation  in  the  values  predicted 
by  these  methods.  The  effect  of  the  diurnal  variations  in  the  photosyn- 
thesis is  seen  to  be  significant  during  daytime  flow  conditions. 
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The  probability  distributions  of  DO  at  different  sections  of 
the  river  shown  in  Figure  20  pertain  to  the  extent  and  intensity  of 
exposure  to  the  daylight.  These  factors  depend  on  the  time  of  initial 
observation  at  the  upper  end  of  the  river  reach  considered. 

In  the  case  of  the  probabilistic  model,  it  is  not  possible  to 
consider  critical  dissolved  oxygen  conditions  as  in  deterministic 
models.  The  objective  in  waste  water  treatment  systems  design  and 
operation  will  then  be  to  minimize  the  probability  of  the  dissolved 
oxygen  level  in  the  stream  falling  below  the  stipulated  value.   This 
could  be  evaluated  knowing  the  probability  distributions  of  the  dis- 
solved oxygen  concentrations  at  different  sections  of  the  river  reaches 
for  the  known  initial  conditions. 

The  writer  became  aware  at  a  late  stage  of  the  possibility 
of  employing  alternative  methods  of  finding  theoretically  the  probability 
distribution  function  for  the  DO  deficit,  when  the  probability  dis- 
tribution functions  for  Kj^  and  K2  and  their  functional  relationship 
with  DO  deficit  are  known.   These  alternatives  are  finding  the  dis- 
tribution of  a  function  of  one  or  more  random  variables  using  change  of 
variable  technique  or  the  moment -generating  function  technique.   The 
analytical  technique  will  yield  a  more  precise  information  than  the 
Monte  Carlo  method. 

A  need  for  taking  into  account  on  a  rational  basis  the  vari- 
ations in  the  parameters  affecting  the  DO-BOD  relationship  had  always 
been  felt  and  this  is  at  least  partly  met  in  this  study.  The  diffi- 
culties one  faces  in  predicting  DO  responses  in  a  river  system  are  best 
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summarized  by  quoting  Kneese  (1964): 

Much  has  been  done  to  develop  generalizations  about 
the  receiving  water  environment  and  to  apply  sci- 
entific principles  to  it.   However  precision  is  less 
and  uncertainty  greater  in  this  area  than  in  most  fields 
of  scientific  and  engineering  forecasting.  This  em- 
phasizes the  importance  of  empirical  checking  of 
forecasts  made  upon  the  basis  of  highly  simplified 
principles. 
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VIII.   ENGINEERING  SIGNIFICANCE 

New  methods  are  reported  for  evaluating  the  composite  effects 
of  combinations  of  waste  water  treatment,  waste  water  flow  regulation, 
stream  flow  regulation,  optimal  allocation  of  stream  dissolved  oxygen, 
etc.  (Liebman  and  Lynn,  1966;  Loucks  and  Lynn,  1966;  Worley,  1963),  for 
application  in  comprehensive  programs  to  improve  and  maintain  the 
quality  of  water  in  major  water  resources  systems.   Several  new  mathe- 
matical techniques  useful  in  fusing  engineering  design,  economic  analy- 
sis and  governmental  planning  have  recently  been  developed  to  exploit 
the  potentiality  of  the  electronic  computer.   In  these  studies,  the 
streams'  self  purification  capacities  and  the  allocation  of  stream  dis- 
solved oxygen  resources  among  various  sources  of  pollution  play  an  im- 
portant role  and  the  reliability  of  the  conclusions  in  these  studies 
depends  on  the  accuracies  of  dissolved  oxygen  predictions,  to  a  smaller 
or  larger  extent.  Thus  the  need  for  estimating  the  response  of  receiving 
streams  to  waste  loads  as  accurately  as  possible  becomes  obvious. 

Kneese  (1964)  advocates  the  need  for  considering  the  proba- 
bilistic character  of  "damage  costs"  which  include  the  optimal  com- 
bmination  of  water  treatment  costs  and  physical  damages  due  to  the  waste 
discharged  upstream  of  the  point  of  use.  The  mathematical  expectation 
of  damages  associated  with  a  particular  level  of  stream  flow  is  obtained 
by  multiplying  the  probability  of  each  flow  by  the  corresponding  cost 
(treatment  plus  damage).   If,  in  addition  to  the  probability  of  occurrence 
of  different  stream  discharges,  the  variations  in  the  basic  velocity 
coefficients  which  define  the  streams*  self  purification  capacities  are 
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taken  into  account,  a  better  estimate  of  the  mathematical  expectation 
of  damages  could  be  obtained*  This  in  turn  aids  in  management  decisions 
aimed  toward  the  best  possible  development  of  a  water  resources  system. 
There  is  a  growing  realisation  for  the  need  of  probabilistic 
oxygen  standards  in  streams  (Loucks,  1965;  Thayer  and  Krutchkoff,  1966). 
In  the  earlier  studies  either  the  chance  variations  in  stream  flow 
(Loucks,  1965),  or  the  variations  in  the  states  of  BOD  and  00  (Thayer 
and  Krutchkoff,  1966)  were  considered.  In  these  cases  the  velocity 
coefficients  K,  and  K2  were  treated  as  constants.  Since  the  variations 
in  Ki  and  K2  have  been  demonstrated  to  be  significant,  it  is  more  realis 
tic  to  consider  the  variations  in  these  coefficients  within  their 
practical  ranges.  The  concepts  developed  in  this  work  can  be  extended 
to  include  variations  in  stream  flows,  thus  obtaining  more  general 
information  on  the  probability  distribution  of  dissolved  oxygen  concen- 
tration in  the  stream.   In  view  of  the  uncertainties  involved  in  the 
stream  self  purification  process,  it  is  more  realistic  to  introduce 
arguments  of  probability  and  predict,  not  how  the  system  will  certainly 
and  invariably  perform,  but  how  it  will  perform  a  certain  fraction  of 
the  time. 
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IX.   SUMMARY,  CONCLUSIONS,  AND  SUGGESTIONS 
FOR  FUTURE  STUDY 


SUMMARY  OF  THE  STUDY 

A  digital  computer  model  is  used  for  defining  the  self 
purification  process  in  the  Ohio  River-Cincinnati  Pool  reach.  The 
model  being  a  representation  of  the  prototype  Incorporating  those 
features  of  the  prototype  deemed  to  be  important  for  the  purpose  at 
hand,  the  operation  of  the  model  in  any  time  Interval  is  character- 
ized by  the  inputs  (BOD)  in  accordance  with  the  parameters  of  the 
process,  yielding  a  sequence  of  outputs  (DO  deficit).   If  the  system 
is  affected  or  perturbed  by  random  components  implicit  in  the  input, 
any  single  simulation  run  yields  only  a  specific  solution  related  to 
one  set  of  conditions  of  the  system.  To  determine  the  general  relation 
of  the  output  to  the  input  and  variables  may  require  a  large  number  of 
simulation  runs.  Monte  Carlo  techniques  are  combined  with  simulation 
analysis  to  remedy  the  lack  of  generality  implicit  in  simulation 
solutions. 

In  applying  Monte  Carlo  simulation  techniques  to  the  Ohio 
River-Cincinnati  Pool  reach,  the  mechanisms  of  BOD  removal  with  con- 
comitant removal  of  oxygen  due  to  bacterial  oxidation  of  organic  matter, 
oxygen  addition  due  to  photosynthesis  and  atmospheric  reaeration  and 
algal  respiration  are  considered.  Other  known  mechanisms  affecting 
the  DO-BOD  relationship  such  as  removal  or  addition  of  BOD  due  to 
deposition  or  scour,  oxygen  demand  due  to  benthal  deposits,  etc.,  are 
considered  insignificant  in  the  oxygen  balance  studies  for  the  particular 
case  investigated.  Oxygen  addition  due  to  photosynthesis  is  treated  as  a 
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periodic  function  having  zero  values  at  sunrise  and  sunset  and  a 
maximum  value  at  mid-day.   Only  arithmetic  average  values  for  algal 
respiration  and  maximum  rate  of  oxygen  addition  due  to  photosynthesis 
are  considered.  The  reaction  coefficients  K^  and  K~  are  treated  as 
random  variables.  A  hypothetical  case  using  Streeter -Phelps  *  formulation 
is  used  to  establish  the  need  for  considering  the  variations  in  Kj  and 

Two  hundred  values  f»r  K.  and  K2  are  generated  by  Monte 
Carlo  methods.   Using  these  generated  values  of  the  reaction  coef- 
ficients and  average  values  for  other  parameters  involved  in  the  process, 
200  values  for  the  dissolved  oxygen  deficit  are  generated,  giving  a 
probable  range  of  values  for  the  dissolved  oxygen  deficit  with  assoc- 
iated probabilities.  The  results  of  the  dissolved  oxygen  deficits 
predicted  by  the  probabilistic  model  are  compared  with  the  observed 
values.  The  predictions  for  dissolved  oxygen  deficits  based  on  the 
most  probable  values  are  found  to  be  in  close  agreement  all  the  time 
with  the  observed  values. 

To  summarize  the  procedure  involved  in  the  use  of  Monte 
Carlo  techniques  for  determining  the  self  purification  capacities  of 
a  stream  under  steady  state  flow  conditions,  the  first  step  is  to  define 
all  the  mechanisms  which  affect  BOD  and  DO,  and  to  formulate  a  mathe- 
matical model  describing  the  DO  response.  Probability  distributions 
for  the  significant  parameters  which  are  found  to  vary  considerably  in 
the  model  are  determined  based  on  actual  observations  for  these  param- 
eters. The  number  of  simulation  runs  required  to  predict  the  DO 
response  under  the  known  initial  conditions  is  dictated  by  the  precision 
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and  the  confidence  with  which  it  is  to  be  predicted.   The  required 
number  of  values  for  the  parameters  which  are  found  to  vary  with 
known  probability  distributions  are  obtained  by  suitable  transformation 
techniques  from  generated  random  numbers.  These  generated  values, 
which  should  conform  to  the  probability  distributions  of  the  historical 
samples,  are  then  used  in  the  simulation  model  to  obtain  probable  range 
of  values  for  the  00  response  with  its  associated  probability 
distribution. 

CONCLUSIONS 

On  the  basis  of  this  study,  the  following  conclusions  can  be 
drawn: 

1.  The  variations  in  the  values  of  deoxygenation  and  re- 
aeration  coefficients  within  their  practical  ranges  have  significant 
effect  in  the  prediction  of  the  state  of  dissolved  oxygen  in  streams 
receiving  organic  waste  loads.   Consequently  there  is  little  justifica- 
tion in  considering  only  the  average  values  for  these  coefficients 
instead  of  taking  into  account  all  the  probable  values  in  any  given 
reach. 

2.  The  error  in  prediction  of  dissolved  oxygen  due  to  the 
variations  in  reaction  coefficients  increases  with  increase  in 
temperature. 

3.  The  values  for  deoxygenation  coefficient  for  river 
samples  obtained  downstream  of  all  the  major  waste  outfalls  in  Cin- 
cinnati, Ohio  are  found  to  vary  at  random.  This  is  probably  due  to 
the  random  changes  in  the  characteristics  of  the  waste  and  the  changes 
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in  the  type  and  distributions  of  microorganisms  responsible  for  the 
stabilization  of  organic  matter.  The  variations  in  deoxygenation  co- 
efficient values  are  adequately  represented  by  Gaussian  distribution 
with  a  mean  value  of  0.173  per  day,  a  standard  deviation  of  0.066  per 
day  and  a  coefficient  of  variation  of  38  percent. 

4,  The  variations  in  percent  error  in  the  prediction  of  re- 
aeration  coefficient  values  for  Tennessee  Valley  rivers  using  regression 
equation  with  mean  depth  and  velocity  of  flow  as  independent  variables 
follow  normal  distribution  law  with  a  mean  of  zero  and  a  standard  devi- 
ation of  36.8  percent. 

5.  The  most  probable  values  for  the  dissolved  oxygen  deficits 
predicted  by  the  probabilistic  model  using  Monte  Carlo  simulation  tech- 
niques are  found  to  be  better  estimates  than  the  values  predicted  by 
the  conventional  deterministic  approaches. 

SUGGESTIONS  FOR  FUTURE  STUDY 

Based  on  the  results  of  this  study  and  the  understanding  of 
the  probabilistic  process  of  the  DO-BOD  relationship  in  streams,  the 
following  suggestions  for  future  research  are  proposed: 

1.  The  applicability  of  Monte  Carlo  techniques  has  to  be 
verified  for  shallow  streams  where  other  mechanisms,  in  addition  to 
the  ones  considered  in  this  study,  affecting  the  DO-BOD  relationship 
are  significant. 

2.  The  interdependence  of  the  reaction  coefficients  K^  and 
Kj  need  further  investigation.  If  they  are  not  found  completely  in- 
dependent, as  have  been  assumed  in  this  study,  a  suitable  joint  density 
function  has  to  be  developed  and  used  in  the  probabilistic  model. 
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3*  The  variations  in  other  parameters  like  oxygen  addition 
due  to  photosynthesis,  algal  respiration,  etc.,  need  to  be  studied 
and  accounted  for  in  the  probabilistic  model, 

4.  The  concept  of  a  simulation  technique  using  Monte  Carlo 
methods  may  be  extended  to  include  variable  initial  conditions  arising 
from  temporal  variations  in  river  flow,  waste  flow  and  waste  strength. 

5.  Research  in  the  area  of  optimal  design  of  waste  treatment 
facilities,  optimal  allocation  of  stream  dissolved  oxygen,  etc.,  should 
include  the  probabilistic  variations  in  the  parameters  affecting 
streams.'  waste  assimilative  capacities. 

6.  The  probabilistic  model  proposed  In  this  study  has  to  be 
solved  analytically  to  determine  the  probability  distribution  function 
for  the  DO  deficit  employing  either  the  change  of  variable  technique 
or  the  moment -generating  function  technique. 
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A}        Coefficients  in  the  linear  regression  equation. 

Ca  Constant  defining  the  critical  value  of  •<!•  statistic  in 
Kolmogorov-Smimov  test  for  samples  of  size  greater  than 
35. 

C  Dissolved  oxygen  concentration,  mg/1,  in  stream  at  any  time. 

Cs  Saturation  DO  concentration,  mg/1,  at  any  given  temperature. 

c.d.f •    Cumulative  distribution  function. 

D  DO  deficit,  mg/1,  at  any  given  temperature  and  time  of  flow. 

Da  Initial  DO  deficit,  mg/1,  in  a  river  reach. 

E^  Oxygen  demand,  mg/(l)(day),  due  to  benthal  deposits. 

D  Critical  DO  deficit,  mg/1,  in  a  stream. 

Dn  Coefficient  of  molecular  diffusion,  L2/T. 

(^(n)     Critical  value  of  'd*  statistic  in  Kolmogorov-Smirnov  test 
for  a  sample  size  of  n  observations  at  a  significance  level 
of  a. 

E(»)  Expected  value  of  a  function. 

ej  Expected  number  of  samples  in  category  i  in  chi -square  test. 

F(x)  Cumulative  frequency  distribution. 

f  Degree  of  freedom  in  chi -square  test. 

f(x)  Probability  density  function. 

g  Number  of  population  parameters  estimated  from  the  sample. 

H  Mean  depth,  feet. 

Kj.       Rate  coefficient,  day"*,  defining  BOD  removal  with  con- 
comitant DO  removal  process  (base  e). 

K2  Rate  coefficient,  day*"*,  defining  atmospheric  reaeration 
process  (base  e). 

K3  Rate  coefficient,  day"1,  defining  addition  or  removal  of 
BOD  process  due  to  scour  or  deposition. 

K5       Rate  coefficient,  day"1-,  defining  BOD  removal  due  to  biological 
extraction  and  accumulation. 
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K_       Rate  coefficient,  day"*1,  defining  overall  removal  rate 
of  BOD  in  a  river  system. 

k  Number  of  observations  less  than  or  equal  to  a  given  value  x. 

L  First  stage  ultimate  BOD,  mg/1,  remaining  to  be  satisfied. 

La  Initial  first  stage  ultimate  BOD,  mg/1,  in  a  river  reach. 

n  Total  number  of  observations  made  in  a  series  of  observations. 

nj  Number  of  samples  observed  in  category  i  of  chi -square  test. 

P{*]     Probability  of  observing  a  value  less  than  or  equal  to  a 
stated  value. 

P|       Number  of  phases  of  length  i  in  "turning  points"  test. 

Pm       Maximum  rate  of  oxygen  addition,  mg/(l)(day),  due  to  photo- 
synthesis. 

Ptt      Rate  of  oxygen  addition,  mg/(l)(day),  at  any  given  time. 

p        Total  number  of  peaks  in  "turning  points"  test  and  period 
of  oscillation  in  a  periodic  sine  function. 

p.d.f.    Probability  density  function. 

R        Average  rate  of  algal  respiration,  sag/ (1) (day),  and  total 
number  of  runs  in  the  "runs  up  and  down"  test. 

Rjj       Number  of  runs  of  length  k  or  more  in  the  "runs  up  and  down" 
test. 

Ro,i>2,.,.,p  Multiple  correlation  coefficient. 

r        Total  number  of  categories  in  chi -square  test. 

z*l       Number  of  runs  of  length  i  in  the  "runs  up  and  down"  test. 

S        Slope  of  river  channel. 

2 
Sg       Sum  of  squared  errors  in  regression  analysis. 

sn(x)     Observed  cumulative  step  function  of  a  sample  of  observations. 

T        River  temperature,  degrees  centigrade. 

t        Time  of  flow,  days. 
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t'       Time  of  flow,  in  days,  of  a  mass  of  water  since  sunrise 
after  that  mass  of  water  has  entered  the  head  end  of  a 
river  reach  considered, 

t_  Flow  through  time,  days,  for  critical  section  in  a  riveri 

V  Mean  velocity  of  flow,  ft/sec,  in  a  river  reach. 

V(«)  Variance  of  a  function. 

X  Random  variable. 

x  Any  particular  value  a  random  variable  X  assumes. 

Zi  Independent  variables  in  regression  analysis. 

ZQ       Dependent  variable  value  predicted  in  linear  regression 
equation. 

Zq'  Actual  value  of  dependent  variable. 

a  Significance  level. 

3  Error  level  in  predicted  values. 

A  Incremental  value  in  the  variables. 

Temperature  coefficient  defining  deoxygenation  rates  at 
different  temperatures. 

$        Temperature  coefficient  defining  reaeration  rates  at 
different  temperatures. 

a       Variance  in  population. 

X2       Chi -square  statistic. 
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APPENDIX  B 

.CTJ  DIAGRAM  FOR  OXYGEN  SAG  COMPUTATIONS; 
SENSITIVITY  ANALYSIS 


€•$ 


Read  temperatures  1  to  30  degrees  centigrade  and  the  corres- 
ponding DO  saturation  value*. 


Read  the  values  of  K,,  K2,  D,,  L,,  river  water  temperature, 
maximum  values  of  K,  (DMK1)  and  K2  (AMK2)  and  Incremental  v< 
In  «!  OBIUKU  and  K2  (AINK2). 


Find  the  DO  saturation  value  corresponding  to  the  river  tempera- 
ture.  Compute  Initial  percentage  DO  saturation  and  write  output 
of  river  water  temperature.  DT.  L«.  Cb  and  Initial  DO. 
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Set  PRESER  -  DEOXK1 
I  -  1 
J  -  1 


De  -  La  exp  (-K^) 


-  r%~  [exp(-K1t)-exp(.K2t) 
K2"*l 

+  Da  exp(-K2t)] 


Call  system 
error 

^--^^  riTN 

NK1 


IK, 


AMK2  .  K2 

ainkV 


+  1.0 


IK,  -  DMK1  '  KW  m 


DINK1 


(K,  -  K.)D 


arc  -  1.0  -  -jqr; 


l'"a 


t  -  0.0 

fi   -  D. 


10      2 
~^~   ln(;r^ARG) 
K2-K1    14 


KlLa 
K2  exp(-K1te) 


D-D  exp(-K1t)+  KjLat  expC-Kjt) 


- 

^C 

-~D~~- 

[0 

+ 

D  «  -s 

DO  -  Cs  -  D 
Percent  saturation 
t  -  t  +  1.0 


C.  -  D 


x  100 


t 

^--^"■~ 

oj 

' 

Write  output  Kj, 
and  critical  DO. 

K2, 

tc, 

Dc 

J  -  J  +  1 


Write  output  sero  to  five  days,  the  corresponding 
deficits,  dissolved  oxygen  concentrations,  and  per- 
cent saturations  for  a  particular  combination  of 
values  for  Kt  and  K2. 
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APPENDIX   C 

FLOW  DIAGRAM  FOR  SOLVING  PROBABILISTIC  MODEL 
FOR   DC-BOD  RELATIONSHIP 
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G> 


(Start) 


Read  input  data;  !)„,  La,  V,  1?,  T,  R,  time  of  initial  observation 
and  other  relevant  parameters. 


Generate  a  sequence  of  uniform  (0,1)  random  numbers  and  evaluate 
the  corresponding  Kj  values.   Correct  the  unreasonably  high  or 
low  values  of  K*. 


Perform  Kolmogorov-Smirnov  test  on  the  generated  values  of  K., 


Is  the  generated  sample  significantly  different  from  the 
historical  sample? 


Yes 


No 


Compute  the  value  of  trend  component  of  K  using  regression 
equation. 


Generate  a  sequence  of  uniform  (0,1)  random  numbers  and  evalu- 
ate the  cor resp aiding  random  components  of  K9  values.   Correct 
the  unreasonably  high  or  low  values  for  the  random  components 
of  K«. 


Perform  Kolmogorov-Smirnov  test  on  the  generated  values  for  the 
random  components  of  K«, 


Is  the  generated  sample  significantly  different  from  the  as- 
sumed distribution? 


Yes 


Develop  random  samples  of  K2  values  by  adding  the  trend  com- 
ponent to  the  random  components  of  K9, 


Set  up  necessary  initial  conditions  for  solving  the  differential 
equations  using  the  sub-routine  for  fourth  order  F.unge-Kutta 
method. 


<D 


Use  a  set  of  generated  values  for  K.  and  K„,  apply  tempera- 
ture corrections  to  them.  Compute  oxygen  deficit  for  the 
known  initial  conditions  using  appropriate  values  for  other 
parameters  in  the  process. 


£ 
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Has  the  last  set  of  generated  values  for  K,  and  K^  been  used?  r.o-h 


i 


Yes 


Go  to  c 


Arrange   the  computed  values   for  the  DO  deficits    in  increasing 
order  and   print    the  results. 


Have  all   input   data  been  processed? 


r;o 


♦  Go   to  B 


Yes 


End 
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APPENDIX  D 

BOD  PROGRESSION  DATA  FOR 
THE  OHIO  RIVER  SAMPLES    (1957  SURVEY) 
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